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SUMMARY: In studies of the stability of aeronautical systems, equations of 
motion are derived which have coefficients dependent on flight speed. Conven- 
tional practice treats the speed as constant, when a set of linear differential 
equations with constant coefficients results. Actually, since the speed varies 
during flight, it may be regarded as a prescribed function of time; the set of 
linear differential equations then has variable coefficients. 

The treatment of the problem of stability then becomes much more 
complex in this case. A simple example is given to show that a system which 
is stable at any constant speed can become unstable during deceleration; the 
ordinary constant-speed criteria are, strictly, therefore inadequate. Some 
approaches to the discussion of stability during acceleration are suggested; a 
solution is given of the single second-order equation which enables the 
amplitude of oscillation of the solution to be studied. Inverse methods of 
approach are suggested, both for single and sets of equations, in which 
particular forms of acceleration corresponding to prescribed solutions are 
derived; and some tentative conclusions are drawn. As would be expected, 
the effects of acceleration depend on a dimensionless ‘‘ acceleration number.” 


Introduction 


The comprehensive study made some 80 years ago by E. J. Routh” has 
provided a firm basis for the assessment of the stability of dynamical systems of 
Which the equations of motion are linear with constant coefficients. Routh’s 

iteria, either in their original form or when modified (as, for example, by Hurwitz 
br Frazer) have been quite remarkably successful in determining the stability of 
actical engineering systems and in pointing to methods of curing undesirable 
stability. It is probably true that most use of Routh’s work has been made by 
heronautical engineers, for whom stability has a peculiar significance: not only 
Hoes the stability of the aircraft itself require study, but provision must also be 
Made against the occurrence of phenomena such as flutter (which may be 
patastrophic) without undue weight penalties. However, many applications of 
outh’s criteria have been made in other fields. 


In aeronautical work, it has always been the practice in the past to treat the 
ir density and the forward speed as constant in stability studies; this has implied, 
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since the coefficients of the variables in the equations of motion are functions of 
speed and density, that these coefficients are constant. Critical (constant) speeds 
for the onset of instability may then be found by the use of Routh’s methods. In 
practice, of course, whether at constant height or not, an aircraft must accelerate 
up to and through such a critical speed; but experience in the past has shown that 
the effect of acceleration, if any, could be ignored. 


However, the very great advances which have been made in the design of 
power plants renders an examination of the effects of acceleration highly desirable. 
In the past, thrusts available have usually been only a small fraction of the aircraft 
weight, and accelerations were correspondingly small (except perhaps transiently in 
diving manoeuvres). Nowadays, however, thrust may considerably exceed weight, 
especially in the field of rocket propulsion, so that accelerations of several “¢” 
may be possible. The question therefore arises: do accelerations of this magnitude 
affect stability, and if so, in what manner and to what extent? 


The present paper aims at drawing attention to this problem, and scratches at 
the surface of some of its facets. Since the speed is regarded as variable with time, 
the problem is basically that of a linear equation or set of equations with variable 
coefficients (the variability may include density or other changes also). However, 
since the stability is in question, it is not necessarily part of the problem to provide 
either exact or approximate solutions to the equations. Such considerations as 
these suggest a variety of approaches to the problem which are not, prima facie, 
obviously closely connected. In setting out the present paper the author has, 
therefore, thought it best to divide it into three Parts. Part I gives some general 
consideration to the problem; Part II deals with some approaches to the solution of 
a single dynamical equation with variable coefficients; and Part III gives one 
possible approach to the study of certain systems governed by sets of such 
equations. 


The various examples given are for the most part in aeronautical idiom, but 
this is purely for convenience and definiteness. 


Part I 


SOME NECESSARY EXTENSIONS OF IDEAS ON STABILITY 


2. What Quantities Determine Stability ? 

As has been remarked already, the equations of motion of an aircraft (or 
perhaps of a component) contain terms dependent on speed, itself a function of the 
independent variable, time; so that strictly the equations are linear with variable 
coefficients. However, since the coefficients usually vary slowly it has been found 
sufficient to treat them as constants, and it has therefore been possible to apply 
Routh’s criteria in the determination of the state of stability over a range of 
(constant) speeds, or, alternatively, to determine a critical (constant) speed at 
which the stability is neutral. 


310 The Aeronautical Quarter) 


Ss 


( 
ti 
fc 
Cr 
t 
d 
ti 
b 
it 
t 
a 
ti 
d 
h 
n 
a 
st 
a 
a 
| | t 
lis 
t 

3 
Ce 
Si 
ti 
t 
| Ne 


STABILITY OF ACCELERATED MOTION 
s of So much work has been done in this way, and the techniques are so well 
ceeds established, that we rarely stop to ask ourselves what we mean by stability in any 
In particular case. In The Principles of the Control and Stability of Aircraft 
rate (Duncan, 1952) stability is defined as follows (p. 113). “A given state of 
that equilibrium at rest or in steady motion of a dynamical system is said to be stable 
when the influence of an imposed disturbing impulse, or force acting throughout 
1 of a finite interval of time, ultimately becomes vanishingly small.” And in develop- 
ble. ment of this theme: “If the solution of these linear equations of disturbance con- 
sate tains no unstable component, the system under investigation is said to be stable 
y in for small disturbances.” This is a familiar conception : we know that with the 
ght conditions postulated all the dependent variables in the equations of motion, and 
Se i their derivatives, become vanishingly small. 
tude But this is not necessarily the whole story; the “influence of an imposed 
disturbing impulse ” may extend beyond the dependent variables and their deriva- 
tives in the equations of motion. For example, the equations of 
. - motion governing the lateral-antisymmetrical motion of an aircraft, unmonitored 
ome by any compass device, do not include the angle of azimuth explicitly; 
able it appears only in its derivatives. Accordingly, if. the aircraft is “stable” 
hen » the derivatives of the angle of azimuth become vanishingly small, but the azimuth 
oe angle itself is arbitrary, as is noted by Duncan (Ref. 2, p. 158). Now if the inten- 
ses tion of the flight is to proceed on some particular heading, then since, following a 
uci | disturbance, the aircraft may settle down on a completely new heading, it can 
has hardly be described as stable; and in order to achieve stability in this sense it is 
eral necessary to use a monitoring device (which might be the pilot) which through its 
n of action on the controls would maintain the heading in flight, and which, from the 
_— mathematical viewpoint, would have the effect of bringing the angle of azimuth 
such J explicitly into the equations of motion. On the other hand, if the criterion of 
Stability is merely “Will the aircraft get into a dangerous or uncontrollable 
_ but condition following the disturbance? ” then the aircraft is in this sense stable. 
This example is given with the object of demonstrating that there is no simple 
_ and universal criterion of stability; we must, in individual cases, consider what 
», are the important attributes of the free flight of the aircraft, and must ensure that 
| these attributes are studied by their explicit introduction into the mathematics of the 
| system. Moreover, it is clear that what is true for equations with constant 
| coefficients will be even more important when the coefficients are variable, since 
| the derivatives of a dependent variable may behave quite differently from the 
variable itself. 
t (or 
fthe | 3. The Wider Problem: Variable Coefficients 
iable | It is a fortuitous but fortunate property of linear equations with constant 
ound |) coefficients that all the dependent variables and their derivatives behave in the 
ipply |} same fashion; and accordingly the state of stability can be studied by an examina- 
e of — tion of the properties of the determinantal equation, without an actual solution of 
d at | the problem. But when the coefficients become functions of the independent 


variable, the first part of the foregoing statement is no longer true. 
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As an example of this, consider a hypothetical aerodynamical system governed 
by the equation 


90c? +3cV x +V?x=0. . (1) 


Here x may be regarded as a disturbance displacement, e.g., normal to a steady 
straight path; c is a typical dimension, e.g., an aerofoil chord; V is the speed. 
When V is constant, the solution of (1) may be written 


=A exp(- cos { B+ 2 


where A, B are arbitrary. At all speeds V, (2) represents a damped oscillation; 
the higher the speed, the greater the damping and frequency. A measure of the 
damping is that, at all speeds, the time to reduce the motion to half amplitude is 
given by 


or Vt= 60c log, 2, 


approximately. Thus the system travels 42 chords while the oscillation drops to 
half amplitude. 


Now, however, let us suppose the system to be decelerated in such a way that 
the forward velocity decreases hyperbolically with time. Then if 


Vt=10c, . ; (4) 
substitution for V in (1) gives 
9x +3tx+10x=0, . ‘ 


of which the general solution may conveniently be written (compare (2)) 


=A (4) cos { (+-)} « « 


where A, B are arbitrary. Thus, when the system is decelerating, it is unstable in 
the sense that x increases indefinitely with time: the derivatives x and x decrease, 
however, as differentiation of (6) will show. In this case a measure of the growth 
of x is that its amplitude is doubled in a time given by 


t 


on the other hand, the amplitude of x falls to a quarter in the same time. During 
this interval the forward distance travelled is 


| Vdt=10c log, + =30c log. 2 
0 
. ‘ ‘ (7) 
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approximately; compare (3). At t=8r, the velocity V has fallen to one-eighth of 
its value at t=t,; the deceleration is thus considerable. The important fact which 
emerges is that a system which, considered as a steady speed problem, is stable 
at all speeds, is shown to be unstable in one particular decelerated motion. It 
follows that we cannot apply steady-speed criteria in the determination of the 
stability of accelerated motions. 


4. Solution of the Equations of Motion ? 


The example given in Section 3 shows clearly that, in contrast to the constant 
coefficients problem, the behaviour of a dependent variable may be quite different 
from that of its derivatives, even though it appears explicitly in the equation. The 
case examined showed that x increased continuously with time while x decreased; 
if the problem is examined in the negative sense of ¢, it is clear that the converse 
is also possible. 


We must therefore ask: do we need to solve the equations of motion to 
provide the time histories of the significant variables and their derivatives? This 
will always give the information required, but the labour involved may be very 
considerable. Methods do, of course, exist for the approximate solution of sets of 
linear equations with variable coefficients: Chapter VII of Elementary Matrices 
(Frazer, Duncan and Collar, 1938) is, for example, devoted to this subject. But 
what is most desirable is a criterion or set of criteria similar to those of Routh, 
by means of which the stability may be assessed without complete description of 
the motions. In view of what has been written here already, it appears that such 
criteria, if found, must determine the behaviour of those variables which are 
important in the phase of stability under consideration and which therefore appear 
explicitly in the equations of motion, and also the behaviour of any important 
derivatives. 


5. Oscillatory Nature of Disturbed Motion 
5.1. Is THERE A CRITERION FOR OSCILLATION? 


There are two aspects of disturbed motion which are of importance in stability 
investigations; they are, of course, related. The more important is whether or not 
the disturbances ultimately die out. But a related matter of importance is 
whether or not the disturbance gives rise to oscillatory motion; some of the most 
important uses of Routh’s criteria relate to the development of unstable oscillations. 


In dealing with an equation with constant coefficients it is a simple matter to 
determine whether the motion is oscillatory or not; it is no longer simple when the 
coefficients are variable, as the following example will show. 


Consider an idealised spring-constrained pendulum, the mass of which is 
effectively concentrated at a distance r from the fulcrum, so that the equation of 
motion is 


which clearly gives rise to a simple harmonic motion. 


November 1957 313 


. 
(3) 
to 
hat 
(4) 
(5) 
in 
ase, 
th 
| 
| 
(7) 


A. R. COLLAR 


Now, however, let r increase with time at constant rate v; then by an appro- 
priate choice of the time origin we have 


r=vt 
This equation has a solution of the form 
6=1? 
where 2p=1+ (1- 
mv? 


Thus, for speeds in the range given by mv? > 4c, the motion is not oscillatory, 
even though the coefficients of 6 and @ in (9) are both positive. For mv? < 40 the 
motion is oscillatory. 


This example poses the question whether we can, even in a simple case, find a 
criterion which will determine whether or not a disturbed motion is oscillatory. 


5.2. SONINE’S THEOREM 


For the very simple case discussed in Section 5.1 a theorem due to Sonine 
(see Szegs™: Orthogonal Polynomials, p. 161) is of help; it is recorded here since, 
if it could be sufficiently generalised, it might provide a starting point for new 
stability studies (some generalisations for a single equation have already been made, 
notably by Polya). 


In its simplest form, Sonine’s theorem is as follows: let 


where ¢ (ft) is positive and has a continuous derivative of constant sign in a range 
R. Consider 


Then, by differentiation of (12) and use of (11) we find 


C= ox’. 2 (13) 
Now the successive maxima of x are defined by x =0, when, by (12), 
|x| 
Thus the successive maxima of |x| will increase or decrease with c. But, 


by (13), ¢ has the same sign as ¢; hence in the range R the successive maxima of 
| x | will increase or decrease with 9. 


Application of the theorem to the example of Section 5.1 is satisfactory: in the 
oscillatory case for t positive, the amplitude increases (in fact, as rf’). But it can 
be seen that the theorem only applies to motion in which x periodically vanishes: 
that is, to oscillatory motion. It is, of course, true that the motion determined by 
(11) is likely to be oscillatory: since ¢(f) is positive, x and x will be of opposite 


314 The Aeronautical Quarterly 


sig 
| 
Mi 
6. 
tri 
os 
wi 

7. 
As 
| 
| 
fo 
WwW 
No 


TO - 


STABILITY OF ACCELERATED MOTION 


signs, so that for x positive the curve will be convex upward and for x negative, 
concave upward: these are characteristic of oscillatory motion (see, for example, 
Mott, Elements of Wave Mechanics, 1952). However, the case discussed in 
Section 5.1 shows that oscillatory motion does not always result from equation (11). 


Part II 


THE SINGLE SECOND-ORDER EQUATION WITH VARIABLE COEFFICIENTS 


6. The Non-Oscillatory Case 


As pointed out in Section 5.1, it is not easy to determine in advance whether a 
solution will be oscillatory or not; probably, in practical cases, it is best to make a 
trial solution, say by step-by-step methods. If, in fact, the solution does not 
oscillate, these methods should be adequate for dealing with the stability of the 
motion. The principal difficulty arises if the solution does oscillate; and attention 
will be restricted to this case. 


7. A Useful Formal Solution 
7.1. EVALUATION OF THE AMPLITUDE FUNCTION 
Let the equation to be studied be written in the form 
Assume as a solution x= AE cos (B+7), (15) 
where €(t), 1 (t) are both real functions of ¢ and A,B are arbitrary constants. 


Substitution in (14) yields the two equations 
Equation (17) is immediately integrable, giving 
=exp (- 2] m dt) 
=exp2, say, ‘ (18) 


The arbitrary constant obtained in the integration of (17) has been taken to be 
unity; no loss of generality is involved, since it is a multiplier of £* and may there- 
fore be absorbed into A. 


From (18), n= | exp 2¢ dt, 20) 


where, similarly, the constant of integration is absorbed into B, 
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Accordingly, a formal solution of (14) is wh 
x= Af cos { B+ exp 20 dr } (21) 
where € is determined by Th 
&+2mé + exp 49, ‘ (22) 
an equation obtained by elimination of 7 between (16) and (18), and where ¢ is oa 
determined by (19). to 
eq 


At first sight, it seems that we have merely replaced the linear equation (14) 
by the more cumbersome non-linear equation (22), of which the left-hand side is — 
identical in form with (14). But if x is an oscillating function, € can be regarded de 
as the amplitude of x and will be a relatively slowly varying function of time; the 


we 
term in cos 7 supplies the oscillations. In these circumstances, (22) will be more 

readily integrable numerically than (14). 

7.2. SOME APPROXIMATIONS, AND CRITICAL CONDITIONS 8. 

In a study of stability, one is normally concerned with those conditions in 
which € is varying slowly, since this is the near-neutral condition. 

Suppose then that (14) gives rise to a slightly damped or slowly growing 
oscillation. Then the derivatives é and £ will be small, and to a first approxima- - 
tion we may neglect them. Then, from (22) - 

and, from (16) or (20), n= : : . (24) 
giving the approximate solution se 

x= An-texp { ~ | mdr) cos {B+ nat } (25) 

It is, of course, possible, by differentiation of (23) and substitution on the left- 
hand side of (22), to obtain a second approximation to €; iterative processes on | 
these lines are readily developed. 

A caution must be added: equation (15) is formally a solution of (14), whether, 
in fact, x oscillates or not. Moreover, any solution of (22) has arbitrary constants, tr 
so that € is neither unique nor always well-behaved. Nevertheless, the solution (15) | 
can be very helpful. 

A critical oscillatory instability is obtained when € is constant, so that the : d 
amplitude of oscillation is neither growing nor decreasing. If we take € to be li 
unity, (16) and (18) give . 

=exp 
so that the critical condition occurs at the time ¢, which satisfies 

n=exp 49, x 
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which, on use of (19), can be written 
AY) n+4mn=0. ; ‘ . (26) 
The use of this expression is illustrated in Section 9. 


2) In the form presented here, the solution (21) was obtained by the author in 
is earlier studies of stability. As might be expected from such an obvious generalisa- 
tion of the solution of the second-order equation with constant coefficients, 
equivalent solutions have been given before, although they do not appear to be 


4) particularly well-known. Hartree“ (Numerical Analysis, p. 148) describes the 
1S F solution for the case m(t)=0, as Madelung’s transformation; Mott (Ref. 5, p. 3) 
~ describes it as the W.K.B. (Wentzel-Kramers-Brillouin) method, but states that it 
e 


was probably first given by Jeffreys. 


re 
8. Some Variants and Linearisations 
in Let x,, x, be any two solutions of (14). Then 
ng 
satisfies ‘y + 6my +2 (4m? + m+2n)y +2 (nr +4mn) y=0, 7 (28) 


as may be readily checked. Now the solution (15) adopted in Section 7.1 
3) is of the form 


4) x, =€ exp (ia), 
so that there will be a second form 

5) x,=£ exp (- ia), 
t- and y=x,x,=& 


will satisfy (28). In fact, the substitution 


transforms equation (22) to 


2yy - y?+4myy + 4ny?=4 exp 49; 


| differentiation of this and elimination of ¢ then gives (28). Equation (28) is a 
e linear equation for y, and hence for &’; its integration by, for example, Frobenius’ 
method may be useful in appropriate circumstances. 


It is to be noted that if €, and hence y, is constant, equation (28) gives the 
condition (26) directly. 


A form of considerable interest can be obtained by the substitution in (14) of 
x=w exp ¢ with, as before, ¢ =—m. 
November 1957 
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The resulting equation is w +[n-m?—m]w=0, 


and has no term in w. All our previous results can be used to interpret this 
equation; in particular, we may write 


where ¢ is determined by the non-linear equation 
. 


Now when the coefficients m, n in (14) are constants, the solution is well known to 
be (compare (2)) 


x=A exp (—mt)cos {B+(n- m?)#t}. 


With m, n as functions of t we now have the parallel form, from (21) on use of (31), 


x= Agexp(-— | ndt) cos { B+ 


where ( is determined by (32); compare also (25). Equation (32) can also be 
linearised by the substitution 
and differentiation. 
While this form has a certain simplicity, the function ¢ is clearly not so useful 
in studies of oscillatory instability as the function é. 


9. An Example 


There is as yet very little experimental evidence on the effect of acceleration 
on stability. In the flutter field, experiments with rocket test vehicles (covering 
accelerations of order 10g and decelerations of order 4g) do not suggest that 
acceleration has much effect on critical flutter speed. On the other hand, it is 
well known that the forward speed at which flutter dies out during deceleration is 
usually appreciably lower than that at which it occurs when the speed is being 
raised; the effect, however, is normally ascribed to friction or to other non- 
linearities in the structural or aerodynamic force-displacement relations. 


The following example is somewhat diffidently offered (it involves some 
assumptions which would be difficult to justify) as an indication of the effects 
of accelerations and as an illustration of the results of Section 7.2. 


The equations of motion of a fluttering system may be typified by 
AX +(BV +D) x +(CV? +E) x+coupling terms=0. . 


Here A is an inertia, BV +D is a damping coefficient of which the first part is 
aerodynamic and the second structural; CV’?+E is a Stiffness coefficient similarly 
compounded. In the absence of coupling terms the system is stable, unless the 
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total stiffness becomes zero, when divergence results. With the coupling terms 
present, however, the system may, at some (constant) critical speed V., develop the 
oscillatory instability known as flutter. 


To discuss this in terms of x only, let us assume that the contributions of the 
coupling terms may be absorbed into the coefficients of x and its derivatives; each 
coupling term will make contributions both in phase and in quadrature, to the 
equation in x, so that the coefficients of x and x will both be quadratic in V. More- 
over, since at the critical speed V. the motion is simple harmonic, the 
coefficient of x must contain a factor V,-V. 


Granted this assumption, we may write the equation in the modified form 


(1 7) a+b 7.) #+ p+q V. x=0, . (34) 
where c is a typical dimension, such as the chord, and a, b,. . are dimensionless. 
Now let us suppose that the system has constant acceleration f, so that 


V=V.+ft, . ‘ ‘ 


the origin of time being chosen at the (constant) critical speed V.. We now 
enquire: at what value of t does the new critical speed occur? 


To find this, let us write 


so that y=1+8, 
and we seek the value of 5. We now apply equation (26), namely, 
n+4mn=0; 


€ 


in this case 2m= — 6 (a+ by) 


pt+qyt+ry 


with i= +2ry 


An important point of general interest to note here is that, for the critical 
speed to change, the stiffness m must be variable. If m is constant, we obtain the 
condition m=0, which gives the constant critical speed. Thus variable stiffness is 
of more significance than variable damping. 


Substitution of the given values of m, n, ri in the condition (26) yields a quartic 
for 5; however, if f is small, 5 will be small, and to first order in 6 we then find 


q+2r+2r) =(), ‘ (36) 
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Then, since K is also small, to first order in 6 and K, 
K (q+2r) 
~ 2(a+b)(p+q+nr)' (38) 


The interpretation of this result is of some interest. We consider the usual 
case, in which the system is stable below (constant) V, and unstable above. In this 
case (a+b) is positive. Since the system is oscillatory, we have also that (p+q+r) 
is positive. Then 6 will have the same sign as K provided that (¢+2r) is positive, 
that is, provided the effective stiffness increases with increasing speed. 


Thus, provided the total effective stiffness increases with increasing speed, the 
effect of positive acceleration is to raise the (constant) critical speed V. by an 
amount V, 5, where 4 is given by (38). 


This result is in accord with expectations based on physical arguments. If the 
effective stiffness grows with increasing speed, and hence with time, then in the 
interchange between kinetic and potential energy, the amplitude of the oscillations 
will tend to decrease; that is, the increasing stiffness has a stabilising effect, which 
is shown by an increase in the critical speed. Conversely, if a fluttering system is 
rapidly decelerated, the energy inherent in the system when oscillating against high 
stiffness may produce very large increases in amplitude if the stiffness is abruptly 
reduced. 


10. An Acceleration Number 


In equation (37) the acceleration appeared in dimensionless form as the 
number fc/V.”._ In general, the simplest acceleration number is probably 


c 


the relation of this quantity to the Froude number is clear. 


11. Inverse Methods 


For dealing with sets of differential equations with constant coefficients, 
inverse methods are often employed. The evaluation of the determinantal equation 
is very laborious if more than three generalised co-ordinates are involved, and very 
considerable reductions in labour are possible by the use of inverse methods. 
Typically, an inverse method is to prescribe a range of solutions which approximate 
to the expected result, and to work back to the corresponding differential equations: 
the optimum match to the actual equations is then chosen. Often the problem is 
reduced in this manner to the solution of a set of algebraic equations for certain 
coefficients which have been left free. 


Clearly, such methods may also be adopted for equations with variable 
coefficients, at least over limited ranges of the independent variable. 
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As an illustration, consider equation (14) and its formal solution (21). The 
function € in (21) is determined by (22), namely, 


£4+2mé +né=£-3 exp 49, 


in which m (and hence ¢) and n are prescribed functions. But it is also possible to 
regard the equation as defining n if m is prescribed and € is chosen arbitrarily; thus 


In this expression, ¢ contains an arbitrary constant, so that for each function € 
chosen a range of curves of n is available for comparison with the actual n. 


To give point to this, consider the case 


m=at—', 
where ais a constant. Then, from (19), 
exp ¢=ct~*, 
where c is an arbitrary constant. If we choose 
where b is also an arbitrary constant, (39) gives n at once as 
c B-—b+2ab 
n= (fa) 


With b,c free for choice in this expression, a wide range of curves of n can be 
drawn for comparison with the curve n prescribed. It may not be possible to 
obtain a good match near t=O but, in view of the expression for m, this region is 
not likely to be critical. 


Part Ill 


ON THE STABILITY OF ACCELERATED SYSTEMS OF AERODYNAMIC TYPE 


12. Sets of Differential Equations with Coefficients Varying in 
Related Ways 


It has been shown already that stability studies of a system governed by only 
one second-order differential equation have considerable difficulty when the 
coefficients are variable. Obviously, when the system is governed by a set of 
differential equations involving coupling terms, the problem is very much more 
difficult: indeed, for the most general type of such a problem, it is probable that 
only numerical methods of solution, giving the complete set of generalised 
co-ordinates and their derivatives, will provide an answer. 
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In the case of certain sets of equations such as arise in problems of aircraft 
stability, however, the coefficients in the equations vary with speed in related 
manners. It is shown in the next section that in such systems, a particular form of 
acceleration enables a general solution to be obtained; the result may be regarded as 
a generalisation of the treatment of sets of equations with constant coefficients. It 
should be emphasised that the method discussed is not of general applicability; it is 
indeed a special approach, by the inverse methods mentioned in Section 11, to a 
set of equations having special properties. Nevertheless, since a complete solution 
is available in this case, it enables certain inferences to be drawn, and brings out 
the relation of stability to stiffness in the manner discussed earlier for a single 
equation. 


13. A Set of Matrix Equations 


Consider the set of equations, in which the order of the matrices is general, 
c?Aq+VcBq+V?Cq=0. . (40) 


Here c is a typical length; V is the airspeed. A, B, C are square matrices of 
inertias, dampings, and aerodynamic stiffnesses; g is a modal column. A set of 
equations such as (40) is the general expression for the equations of motion of a 
rigid body subjected to air forces. With the equations in the form (40), A, B,C 
and q can be dimensionless. 


When V is a constant, say v, we solve (40) by writing 


Gq =(*)'q. . (41) 


Then substitution in (40) gives 
. ‘ ‘ (42) 


and the permissible values of A are given by the determinantal equation 


When V is not constant such a solution is no longer true. However, let us 
suppose we can find a solution of (40) with V variable in the general form 


where s is a scalar function and q, is a column of constants. Substitution in (40) 
then gives 
['s c?A+ VcB+sV°C] q,=0, 
and the elements in the square matrix are functions of ¢. If, however, we 


assume Sc? and sVc to be proportional to sV?, this function will cancel, and we 
shall have a square matrix of constants. Accordingly, let us write 


Sc=msV; gc=nsV, 
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where m,n are constants: they are not entirely arbitrary, since we have a deter- 
minantal condition to satisfy. Elimination of V/c between the two equations 
enables us to write 


ae 
which integrates to s=(a+ bt)’, 


where a, b are arbitrary and r is related to m,n. Finally, V is inversely propor- 
tional to a+ bt. 


Accordingly, we return to equation (40) and write 


so that at t=0 we have 
V=v, 
dt 


In view of the derivation of the hyperbolic form for V, we shall therefore have a 
solution of (40) of the form 


Substitution of (45) and (47) in (40) gives, after cancellation of the common function 
of time, 
[c’?f?r (r- 1) A—v°cfrB + v*C] q,=0. 


where K is an acceleration number, as given in Section 10; then 


[A (A+ K) A+AB+C] q,=0 


and A is given by the vanishing of the determinant of the square matrix in (49). 


14. Comparison with the Steady-Speed Case 
If in equation (45) f is made vanishingly small, V becomes constant at the 
value v. Equation (49) reduces at once to (42). Moreover, on substitution from 
(48) for r in (47), we obtain 
_ ft —A/K 
a= (1-5) 


which, in view of (48), may be written 
it -A/K 


November 1957 323 


aft 
ed 
of 
as 
It 
is 
on 
le 
al, 
0) 
of 
of 
C 
) 
1S 
‘ 
4 
| 


A. R. COLLAR 


and, when K becomes vanishingly small, this expression reduces to 


q=o 


compare (41). Thus the solution (47), corresponding to the hyperbolic variation 
of V with time given by (45), is a generalisation of the usual steady-speed solution. 


15. Values of A for Stability 


We assume ¢ always increasing from ¢=0, that is, t positive; also v is positive 
for systems of the type (40). 


Consider first a positive acceleration, that is, K positive. Then the bracketed 
term in (50) reduces linearly with time from 1 to 0 as V rises from v to infinity. 
If, therefore, q is to reduce in magnitude, the index in (50), or its real part if it is 
complex, must be positive; this implies 

KR (aA) <0, 
where & (A) is the real part of A. 


Conversely, when K is negative, the bracket is positive and exceeds unity. 
Thus the index in (50), or its real part, must be negative, and again we find 


RA) <0. 


16. Evaluation of the Stability 

Write equation (49) in the form 

The equation is now in the standard form; and since for stability we require only 
the standard result 
RA) <0, 

we may treat (51) by means of Routh’s stability criteria. 

We thus find the simple rule that, to investigate the stability of the system (40) 
when V is varying according to (45), we add to the damping matrix B the inertia 
matrix A multiplied by the acceleration number K; the system is then used to 


derive the determinantal equation implied by (51), which is then examined by 
Routh’s criteria in the usual way. 


17. An Example 


Consider the system 


Bla 106, 1007. 7186, 1807 _ 

40, 38 70, | . (52) 
v? 

with V= (53) 
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Case 1. No Acceleration 


With f=0, we have V=v and 


which leads to 


The determinantal equation becomes 

A* + 164° +54A7+761+48=0. . (56) 


Application of Routh’s criteria shows stability: for the coefficients in (56) are all 
positive and 


T,=16 x 54 x 76 — (76)? — (16)? x 48 = 47,600 > 0. 
In fact, (56) may be factorised to 
(A.+ 12) (A+2) (A2+2A4+2)=0, 


whence it follows that the complete solution of (52) with constant V =v is 


(57) 
where k,, k,, k, and ¢ are arbitrary. 
Thus at constant V the system is stable. 
Case 2. Deceleration 
Consider the case for which K = - 4; by (45) this gives 
V c 


Let us examine this case by an application of Routh’s criteria. We first 
evaluate 


106, 100 5, 27 786, 92 
and (51) then gives 
5A? + 86A+ 186, 2A?+92A+ 180 
10, A? +34A+ 6g, | + 69 
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The determinantal equation is now 
At + + 6A? — 28A + 48=0, . (60) 
and Routh’s criteria show at once that the system is unstable. 
Again we may factorise the determinantal equation: it is 


(A+ 6) (A+ 4) (A?-—2A4+2)=0. 


We are now able to evaluate the complete solution of the system of equations (52) 
with V given by (58), namely 


It is 


where k,, k,, k, and < are arbitrary. The last (oscillatory) term in this expression 
grows indefinitely with time, confirming the result found from Routh’s criteria. 


18. The Critical Acceleration for Instability 
The example of Section 17 illustrates the fact that a system which is stable at 
a steady speed may become unstable under conditions of acceleration. 


Routh’s criteria can, of course, be used in this form of analysis to determine 
critical conditions. Just as, in a flutter problem, the criteria can be used to find 
critical speeds, so in this problem we can find critical accelerations; for the 
determinantal equation implied by (51) is a polynomial in A with coefficients which 
are functions of K, so that use of Routh’s criteria can provide the critical values 
of K. 


Let us illustrate this by the example of Section 17. With K general, 
equation (51) yields the determinant (with the matrices given by (52)) 
+(106.+5K)A+ 186, 2A?+(100+ 2K) A+ 180 
207+ (40+2K)A+ 70, (38+ K)A+ 


or A* +(16+ 2K) AS + (54+ 16K + A? + (76 + 26K) A+48=0. (62) 


The conditions for stability, in this case, are therefore that the three functions of K 
in (62) shall be positive, and that T, shall be positive, ;where 


T,=(16 +2K) (54+ 16K + K*) (76+ 26K) —(76+ 26K)? - (16+2K)? 48 
= 52K‘ + 1400K* + 12244K? + 43104K + 47600. . (63) 
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Now T, has the factors (K +2), (K.+14), (13K? + 142K + 425), and accordingly the 
first instability of the system, which is clearly stable for all positive values of K, 


| occurs when 7, vanishes for K=—2. Hence, for 


f=-Q'/e, . . (64) 
the system will be neutrally stable. 


19. Lateral Stability of an Accelerating Aircraft 


The system discussed earlier is special in character, since all the stiffnesses 
are of aerodynamic origin. The introduction of elastic or gravitational stiffnesses 
would require another approach. 


The standard equations of lateral-antisymmetric motion of an aircraft in 
horizontal flight are 


m 
Ap pL,—rL,=0 
pN,-rN,=0 

with =p. 


' If we wish to discuss conditions in which V is variable we cannot make all the 


usual substitutions to derive dimensionless quantities, since V must be excluded as 
a factor. Accordingly, we write 


 sii=v; 


here s is a linear dimension. The equations can then be written in the form 


0, -iz, ic| | 0, 0,0} 

. 


where i4, Yy, #2, and so on, have their usual meanings. Finally, if we substitute 
gt?=sT’ 
gs, 
so that T, N are dimensionless measures of ft, V respectively, we obtain 
Aq’ +NBq’'+Cq=0 . ‘ (66) 
where A, B,C are the square matrices in (65), 


and primes (’) denote differentiation with respect to T. 
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It is to be noted that, strictly, the coefficients i,, ic, ig and certain of the 
derivative coefficients are functions of C, and hence of N; however, we shall 
assume that any variation in N is such that variations in the coefficients can be 
neglected, and we shall treat A and B as constants. 


As before, we look for a solution of the form 


q=2(t)q, 


where 2 is a scalar function and q, is a column of constants; and we see that such a 
solution is possible if 
a”, No’, a 


are all proportional to the same function. This condition is satisfied by, for 
example 


a=sinh AT 
(67) 
N=n tanh AT 
and then we require 


so that A determines the exponential roots corresponding to a fixed speed N=n. 


The interpretation of the solution 


sinh AT 
N=ntanh AT 


(69) 


is of interest. Thus:— 


(a) Let A, be positive; that is, it is an unstable root of the system with constant N. 
Then (69) shows that g and N grow with T from T=0. When T is large, 
N tends to the constant value n and g to the usual exp(A\7) form. Thus, with 
the system accelerating, it is still unstable. 


(b) Let A, be negative. For physical reality we require both N and n to be 
positive; accordingly we consider only negative values of T (the origin of time 
is, of course, arbitrary). We then find that N is decreasing from the constant 
value N=n towards zero; and similarly q tends toward zero. Thus a stable 
root remains stable under decelerating conditions. 


(c) The other constituent motions do not appear to be simply derivable when N is 
varying in the particular way (69) corresponding to the root A,: the accelera- 
tion can be made to correspond to any one root at a time and the related 
constituent is then given by (69), but it does not appear easy to discuss the 
effects on the remaining constituents. 
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(dq) An alternative solution is 


(70) 


q=q, cosh 
N=ncothaT 


In this case, if A; and T are positive, the system is decelerating, but is still 
unstable. 


(e) If A, is negative we consider the negative range of T: then as T increases 
toward zero N becomes indefinitely large, but q decreases; that is, the system 
remains stable. 


Thus we find from (a) and (d) that if there is an unstable constituent A, at the 
fixed speed N=n, it remains unstable whether the speed increases toward n (as 
tanh A; 7) or decreases toward n (as cothA, 7). Similarly, from (b) and (e), a stable 
constituent remains stable under similar acceleration or deceleration. We cannot, 
however, deduce the effects of these changes of speed on the other constituent 
motions. 


20. Some Comments 


The system of equations (65), being a set of three second-order differential 
equations, has a determinantal equation with six roots. The matrix C is, however, 
doubly degenerate, so that two roots are zero. One of these (see Ref. 2) denotes 
neutral equilibrium in azimuth; ¥ is not determinable. The other denotes neutral 
equilibrium in the quantity », the measure of lateral displacement. 


In some respects the treatment of the system of Section 19 is less satisfactory 
than that of Section 13. It was possible for the first system to apply Routh’s 
criteria to determine critical conditions and accelerations; only a relatively limited 
examination of the system of Section 19 was possible. This may stem, in part at 
least, from the choice of axes; and it may be that the “wind axes” fixed in the 
aircraft, which have provided the framework on which have been built studies of 
stability at constant forward speed, will prove to be less suitable for studies 
involving acceleration effects. 


21. Conclusion 


It is clear from the foregoing paragraphs that the study of the stability of a 
dynamical system, for which the differential equations have variable coefficients, 
is much more difficult than with constant coefficients. To begin with, some 
extensions of the concepts of stability are involved; the stability is no longer 
governed by exponential modes to which a variable and all its derivatives conform. 
Criteria involving all relevant variables and all their important derivatives must 
be sought. 


At the same time, complete solution of the equations is evidently very 
laborious in general and may often be impossible except by numerical methods. 


November 1957 329 


| 

67) 

N. 
ge, 
ith 
me 
ant 
ble 
] is 
ra- 
ted 
terly 


A. R. COLLAR 


It is possible to envisage cases in which a solution must be examined over a very 
wide range of the independent variable because of an apparent incipient instability 
which might in the event never develop. Criteria for stability similar to those of 
Routh are therefore of great importance. 


It seems from the studies given here that inverse methods may be important 
in approaching this problem. The inverse methods illustrated in Part III 
undoubtedly throw light on tendencies and direct attention to important parameters. 


As would be expected, acceleration effects depend on a dimensionless 
parameter which appears to arise naturally in the analysis in the form of an 
“acceleration number” K=fc/V*. In practical cases it seems probable that K 
will not be large; if very high thrusts are available to give high f, correspondingly 
high values of V will also be probable, except near the beginning and ending of a 
flight. For small K, it seems probable that acceleration effects will not be very 
significant (see, for example, Section 13). 


Finally, the studies both of Part II and Part III appear to indicate that, for 
acceleration effects to be important, the stiffnesses appearing in a problem must be 
variable. In Part II it was found that variable damping alone did not alter a 
critical speed; it was necessary that the derivative of the stiffness should not vanish. 
In Part III, a problem involving (variable) aerodynamic stiffnesses gave significant 
acceleration effects; a problem with (constant) gravitational stiffness showed no 
stability change with acceleration. Such changes as occur seem to suggest that it 
is the interplay of kinetic and potential energy that determines the stability; where 
the stiffnesses are varying the potential energy for a given amplitude also varies. 
If this is true, it may be that a search for a criterion based on the “energy of 
disturbance ” will be most profitable in stability studies of this kind. 
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very} Introduction 


This paper gives a brief account of the more detailed experiments which were 
t, for used for studying the “ Jet Flap” mechanism. Previous papers“: * » and a lecture“? 
have dealt with the theoretical and other aspects of the subject and it has been 
shown that, for an uncambered aerofoil at zero incidence having a thin full span jet 
issuing downwards from the trailing edge at angle @ to the chord line, the lift 
coefficient is substantially 
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- where C,;= jet coefficient 
aries. 

gross thrust (J) 


~ x (wing area) 
Another useful concept is the lift gain, 


nica- Ga Total lift 
Nos. = Tet lift 


which, in its simplest form, is given by 


ersity G=5/C,}. 
| Following an early crude experiment” which roughly confirmed the lift hypothesis, 
ridge § it was important to prove the thrust hypothesis also, as well as to provide detailed 
information on the lift and pitching moment both with and without incidence. In 
addition, more detailed investigations were made as required by the growing under- 
standing and some of these also are described in this paper. To save space, 
| illustrations already used in other papers of the series"’<* are referred to but not 
reprinted; a more detailed account and all the experimental results appear in 

| Reference 5. 
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NOTATION 
U, undisturbed velocity relative to aerofoil 
U, velocity in the jet 
Py Main stream density (assuming incompressible flow) 
ps density of jet fluid 
R_ chordal Reynolds number 
aerofoil chord 
a angle of incidence 
n analogous flap angle 
jet deflection angle relative to chord line 
J total jet reaction or momentum flux at the nozzle 
Cy, “no blow” drag coefficient, i.e. when C;=0 
Cy ets effective drag coefficient =C, - Cr, 
Cy, thrust coefficient, aerofoil at zero incidence 
Cp, jet drag coefficient 
C, jet coefficient 
C,, total lift coefficient, aerofoil at zero incidence 
C, pressure coefficient 
G lift gain factor 


2. The Apparatus 


The model aerofoil was limited in size by the air supply available, but the 
chord and thickness were required to be as large as possible from mechanical 
considerations. Therefore, a chord of 8 in. was chosen with a span of 12 in. and 
a thickness/chord ratio of 124 per cent. The ideal profile for a jet flap wing was 
unknown, and so an ellipse was thought satisfactorily neutral, with the advantage 
of being simply transformable to a circle should the trend of theoretical work require 
it. An external view of the model appears in Ref. 4, Fig. 4. This size of model 
permitted 26 static holes to be fitted and a section of the aerofoil with static hole 
positions is shown in Ref. 4, Fig. 6. In addition eight static holes were added at 
the mid-chord station and distributed across the span to check that effectively two- 
dimensional flow was obtained at high lift conditions. 


Already it had been decided that, in order to obtain a direct answer to the thrust 
problem, the first model aerofoil should have its full span jet slot deflected 90° from 
the chord line and that alternative trailing edge assemblies should be fitted for other 
deflection angles (nominally 30° and 60° in order of testing, the nominal angle 
being used here to designate the model). 


A wind tunnel width of 12 in. being determined by the model span, a deep 
working section was chosen to avoid the interference effects probable with the 
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expected large lift coefficients. Therefore, a two-dimensional contraction maintained 
a depth of 564 in. from the square outlet of an existing low speed cascade tunnel and 
the arrangement is shown in Ref. 4, Fig. 3. An air speed of 100 ft./sec. was attain- 
able, with a small margin in reserve, and the velocity distribution across the working 
section was found to be reasonably even. 


The air for the jet (about 0-17 Ib./sec. maximum) was supplied through the 
100 Ib. /in.? laboratory ring main by a diesel driven compressor. Two-stage pressure 
reduction was used to minimise the fluctuations due to the compressor governing 
mechanism, but the manual control of a blow-off valve was still required at small 
air flows (small values of C, in general) and also for large air flows beyond the 
control capacity of the reducing valves. It was found important to remove con- 
densed water from the air, as any passing through the jet slot could give rise to large 
erratic variations in thrust and a noticeable effect on the pressure distribution. A 
continuous “ blow-down ” of condensate and oil from all the drains in the lines and 
separators prevented this trouble and the amount of liquid removed amounted, on 
occasions, to more than two gallons per hour. 


The air was fed into the model at both ends through a tubular “ main spar ” from 
vertical lengths of rubber tubing which were always adjusted to be straight and taut 
to avoid any Bourdon tube effect when under pressure during the thrust measure- 
ments. Small holes along the “spar” throttled the air into the body of the 
aerofoil, so promoting an even spanwise distribution, while the area contraction 
from this point to the jet slot was more than 40 to 1. 


A thrust balance was constructed with a simple parallelogram linkage on either 
side of the tunnel, the forward pair of arms carrying a torque tube with the weigh- 
beam fixed to it centrally, and a sensitive null-point pneumatic indicator made the 
measurement of thrust and drag forces reliable to within 0-002 lb., while a change 
of 0:0005 could be detected. 


Some difficulty was found in ridding the model of leaks which, if not prevented, 
would have upset its boundary layers. The strength of the model limited the jet 
total pressure (measured arbitrarily within the body of the aerofoil) to 15 Ib./in.? 
corresponding to a C, of about 0°5 at U,=100 ft. /sec. 


Figure 1 (p. 341) shows the model in the working section. 


3. The Basic Experiments 


3.1. MODEL CALIBRATIONS 

With the chord line marked on the end of the model and a horizontal line 
scribed on the transparent tunnel wall, checks were made that the air flow in the 
empty tunnel was in fact horizontal and that zero geometric incidence of the model 
was also the angle of zero lift with no jet. 


An accurate assessment of the angle # between the mean jet centre line and the 
chord line was essential, the degree of accuracy required increasing as @ decreased. 
This was because the pressure thrust (or drag) could only be measured as the 
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The lift at zero incidence. 


difference between the directly measured force and the (calibrated thrust) x cos 6, 
unless an integration of the static pressures was to be relied upon. After trying an 
attractive null-point method, which failed because of the resultant effects of the jet 
impinging on the top wall of the wind tunnel, the most reliable method was found 
to be the simplest: that of observing fairly long, fine wool strands attached to the 
trailing edge, close to and on either side of the jet slot. 


For the thrust calibration, the model was set at an incidence of -9@ to give a 
horizontal jet and the balance readings then obtained over the full range 
(0-15 Ib./in.”) of jet total pressure were corrected by adding on the pressure drag 
caused by the entrainment flow around the model. This “ induced ” force naturally 
varied with the attitude of the three models and the corrections amounted to about 
2:5 per cent, 6:1 per cent and 9 per cent for the 30°, 60° and 90° models respectively. 


3.2. STANDARD TESTS 


For all three models the minimum information required was the measurement 
of total lift and pitching moment, both at zero incidence and at incidence, and the 
resultant thrust (or drag) by balance at zero incidence only. The thrust balance was 
not complete in time for the 90° model tests but a “ jury ” rig gave some rough and 
ready measurements for comparison with those obtained by pressure plotting. 


The results are given in Ref. 5, the measured total lift coefficients for the three 
basic models at zero incidence being correlated here in Fig. 2 in the form C,,/sin 4 
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Aerodynamic centre and centre of lift. 


plotted against C,;!. The theoretical curve from equations (6) and (10) of Ref. 1 
is also added, as is the result of the simplified formula C,,=5 C,' sin 0. 


The centre of lift position (zero incidence) for the whole range of C,; can be 
seen in Fig. 3, where the theoretical curve is that for »/sin 6=1 in equation (13) of 
Ref. 1. However, even had the approximation been »=@, the changes in the 
theoretical curve would have been insignificant beside the discrepancy from the 
experimental points. Figure 4 illustrates the thrust measured at zero incidence, both 
by balance and by pressure plotting, with the results of the 0° model tests added 
(see Section 4.4). The existence of a pressure thrust is undeniable but its magnitude 
relative to the full jet reaction was initially rather disappointing. However, on 
investigation, so many sources of “loss ” became apparent that the pessimistic view 
was not maintained (see Section 4). Examples of the pressure distribution for the 
three models are shown in Fig. 5. They are the distributions corresponding to a 
total lift coefficient of approximately 2:5 for all the models, the change in 6 and C, 
producing the expected modification to the peak suctions at the leading and trailing 
edges. The results from the tests with the aerofoil at incidence are illustrated only 
in Fig. 6, although curves of C, against the incidence, 2, will be found for the 30° 
and 60° models in Figs. 17 and 19 of Ref. 4. 
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The measured thrust with a deflected jet. 


4. Other Investigations 


4.1. SOME REYNOLDS NUMBER EFFECTS 


For all the “standard ” tests with a jet coefficient of 0-5 or less the chordal 
Reynolds number, R, was 4-25 x 10° (main stream speed = 100 ft./sec.) but, with the 
pressure inside the model limited to about 15 Ib./in.? gauge, values of C, above 0°5 
could only be obtained by reducing the air speed in the tunnel, and so also the 
Reynolds number. The tests on the first model, with a jet deflection of 90°, 
reached an upper limit of C,;=0-5 since, at that time, it was thought unlikely that 
such a value would be exceeded in practice. 


Analysis of the 90° model tests prompted some experiments with transition 
wires fitted to the aerofoil, which, by then, had been rebuilt with its 30° trailing 
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edge. As a result of all the observations made, the following interpretation is 
suggested for these experimental conditions and illustrated in Fig. 7. 


With values of C, sufficiently small to avoid the pressure coefficient near the 
leading edge exceeding a certain critical value (about - 1-0 for these models), a 
laminar boundary layer persists nearly to the trailing edge when no transition wire 
is fitted and then leaves the surface to give a large loss. This effect is indicated by 
the lack of pressure recovery on the upper surface near the trailing edge of the 30° 
(and 60°) model (the full line of Fig. 7(a)) and a correspondingly low lift (lift is 
plotted as lift gain G to emphasise the effect). With transition wires fitted there is 
an increase of static pressure around the trailing edge (the dotted line of Fig. 7(a)) 
and with the upper surface wire near the trailing edge at “C,” the greatest benefit is 
obtained from having the laminar boundary layer from the leading edge to the wire. 
When the peak C, near the leading edge exceeds - 1:0 the subsequent turbulent 
boundary layer does not separate at the trailing edge, except for the 90° model, 
where separation occurred at all times after a partial recovery following the 
forward transition. 


At the lower chordal Reynolds number of 1:28 x 10° (flow visualisation tests, 
60° model) the peak negative pressure coefficient for transition near the leading edge 
appeared to be above 2:4. 


Leading edge separation occurred with the 30° and 90° model after a peak 
pressure coefficient of about - 6-7 had been reached (C,=1-50, R=2°12 x 10’, 
30° model; C,;=0:38, R=4-25 x 10°, 90° model) and for the 60° model the peak 
suction was C,= -— 7:5 (C,;=0-50, R=4-25 x 10°) before the flow separated, but this 
higher value might have been due to the thinner boundary layer on the tunnel walls 
obtained after improvements had been made to the joint between the contraction and 
the working section. An example of the pressure distribution before and after the 
onset of leading edge separation is given in Fig. 8. 


4.2. FLOW VISUALISATION 


A simple smoke generator was constructed, similar to that developed by Preston 
and Sweeting“ in which kerosine was boiled and the jet of vapour condensed by 
two cold air jets to produce a dense mist. The “smoke” was injected into the 
tunnel in ten filaments and illuminated from the open tunnel exit. 


At a tunnel speed of 30 ft./sec. reasonable delineation of the smoke trails was 
found over the full range of C; from 0 to 5-0. With no transition wires fitted to the 
model and with the low Reynolds number of 1:28 x 10°, it was obvious that the 
recorded flow patterns could not be compared directly with the performance 
measured at the higher Reynolds numbers of the “ standard” test. For this reason 
pressure distributions and thrusts were measured at the same Reynolds number and 
at all values of the jet coefficient C,; at which photographs had been taken. It was 
found that the pressure lift was about 7} per cent less than for the tests at higher 
Reynolds numbers, that transition near the leading edge occurred at a lower value of 
C, but a larger (negative) value of C, and that leading edge separation appeared at a 
lower value of C;. The photographs used to illustrate Ref. 4 are sufficient in 
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Fig. 9. Crown Copyright 
A typical flow pattern for the 60° model. 


40 


10:75 


1-28x108 
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The flow in the region of the trailing edge. 
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Fig. 10. 


Ground interference effect on lift with the 30° and 60° models. 


number to obviate an addition here, but a more general view of the smoke technique 
is illustrated in Fig. 9, the photograph being taken from the other side of the wind 
tunnel after a large transparent panel had been fitted for cinematography. It should 
be noted that, with “no blow,” the second smoke filament from the top would be 
divided by the model. 


4.3. GROUND INTERFERENCE 


Pressure plotting tests at zero incidence only were performed on the 30° and 
60° models and thrust measurements were confined to the three smallest ground 
clearances and the 60° model. With the 30° model, the effect of ground proximity 
was far from excessive. The reduction in lift and thrust with the 60° model, how- 
ever, was found to be rather large, particularly at a ground to chord line distance of 
only 3; of the chord. Figure 10 illustrates the effect of ground interference on lift, 
while the full results will be found in Ref. 5. 


The three-dimensional effects must have been large with the small aspect ratio 
plus wide end gaps caused by the “ ground ” being wider than the model span. The 
induced drag in particular must have been considerable for the 60° model with its 
large lift coefficients (nearly twice those of the 30° model). The results of the 
ground interference tests may therefore be pessimistic, but repeat tests with an 
improved “ ground ” could not be undertaken in the available time. 


4.4. DRAG WITH AN UNDEFLECTED JET 

While examining the air flow around the 90° model with a smoke probe and 
wool tuft, it was observed that the main stream was entering the jet sheet near the 
aerofoil trailing edge at very large entrainment angles—those at the lower edge of 
the jet being more than 90° to the jet centre line over the first inch or so. Figure 11 
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The sink drag effect for the 0° model. 


shows this effect with the 60° model. These large entrainment angles suggested a 
“loss” in the mixing process and prompted the construction of a simple aerofoil 
having an undeflected jet slot situated in a comparatively thin trailing edge 
(Fig. 24 of Ref. 4). 


An “effective ” drag coefficient, Cp o:, was found, using a cold air jet, by sub- 
tracting the coefficient of measured thrust (C;,) from the jet coefficient (C;). The 
result is shown in Fig. 12, with the results of two additional tests; the one made with 
a spoiler fitted to the model leading edge to increase artificially the “ no-blow 
drag coefficient, Cp, and the other with an improved leading edge shape of 
“ Plasticine ” which gave a lower Cp,. It can be seen that the slope of all these 
Cp et: -— C; curves is about the same, thus indicating that the jet drag, or sink effect, 
is a pressure drag caused by a low pressure region in the vicinity of the aerofoil 
trailing edge, as expected, and that it amounts to about 6 per cent of the gross thrust. 
By subtracting Cp,, all these data can be reduced to the common basis of a jet drag 
coefficient, Cp;, which has only a rough practical value, since it includes the extra 
skin friction due to local increased velocities. 


” 


Figure 9 gives an indication of the contraction of the main stream which occurs 
when p,U; > p,U,*”, even when due allowance is made for the effect on the 
streamlines of the potential flow distribution. 


The 0° model was also used with a hydrogen jet to investigate the effect of the 
important factor p,U,/p,;U,; for a jet-flapped aerofoil under cruising conditions. 
Considerable difficulties were experienced in the experiments, but the main result 
supports Stratford’s theory and is quoted and discussed fully in Ref. 2. 
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4.5. THRusT Boost 


The mixing theory of Refs. 2 and 3 gave rise to the prediction of thrust 
augmentation under ideal conditions and an experimental confirmation of this effect 
was therefore attempted. Unlike the 0° model, however, the thrust boost model 
was quite refined and very well made, as can be seen in Fig. 4 of Ref. 2 (also Fig. 28 
of Ref. 4). The “fish tail” exhausts for the main jet were individually calibrated 
and adjusted to pass the same flow for a given reservoir pressure to avoid major 
secondary flows across the span, with consequent drag. An accurate slot, 0-005 in. 
wide, was provided at the hinge of the flap and separately supplied with air for 
boundary layer control, should the mixture of jet and main stream prove recalcitrant 
at this point. In practice, however, this was found to be unnecessary and the 
main jet sheet negotiated the “corner” without any overall separation, although 
see Ref. 2. The experiment is described and discussed fully in that paper (and 
briefly in Ref. 4) and so needs no repetition here. 
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Noise Measurements on Cold Jets using 
Convergent-Divergent Nozzles 


R. E. FRANKLIN, Ph.D. 


(Department of Aeronautical Engineering, University of Southampton) 


SUMMARY: Tests have been carried out on a series of convergent-divergent 

nozzles and the results compared with those from a test on a convergent nozzle. 

At the design pressure ratio, reductions in the noise were found, particularly 

at 90° to the jet direction. For lower forward pressures the noise is also less, 

but above the design pressure the noise rises rapidly to the same level as that 

with a convergent nozzle. In the experiments cold jets were used, but the 
results are probably applicable in certain respects to hot jets. 


1. Introduction 
1.1. FLOW CHARACTERISTICS 


The behaviour of a gas as it flows through a convergent nozzle is best under- 
stood by first considering the case when the gas flows from a reservoir, in which the 
conditions are constant, through the nozzle into a large receiver in which the static 
pressure can be varied. At first, as the ratio (pressure in reservoir)/(pressure in 
receiver) is increased from unity, the mass flow through the nozzle increases, the gas 
issuing into the reservoir in a uniform stream. But the increase in mass flow 
corresponding to a given increment in the pressure ratio gradually decreases until 
the pressure ratio reaches a critical value at which the rate of change of mass flow 
becomes zero, and further increase of the pressure ratio will not affect the mass flow. 
The explanation of this is to be found in the fact that at the critical pressure ratio 
the velocity of the gas in the exit plane of the nozzle becomes equal to the local 
speed of sound, and consequently further reductions in the receiver pressure cannot 
be felt inside the nozzle. The pressure ratio of the nozzle-controlled expansion 
thus stays the same for all overall pressure ratios greater than the critical. It is 
usual to say that in this condition the nozzle is “ choked ”. 


Although for pressure ratios above the critical the mass flow through the 
nozzle does not change, the character of the issuing jet undergoes a great change. 
When the nozzle is choked the gas crosses the exit plane of the nozzle at a pressure 
greater than the receiver pressure. Further expansion must take place and this 
begins at the nozzle lip, from which springs a fan of expansion wavelets which 
penetrates into the jet, causing the jet to diverge. These expansion wavelets are 
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NOISE MEASUREMENTS 


reflected at the constant pressure boundary as compression waves and the jet 
contracts again so that at a certain distance from the nozzle the conditions in the 
exit plane are repeated. Theoretically this is repeated indefinitely and it is seen 
that the flow in the jet is cyclic. 


This type of flow has several disadvantages, some of which will appear in what 
follows, and may be avoided by fully expanding the gas down to the receiver 
pressure. To do this, the gas must be further accelerated to supersonic speeds and 
this may be done by adding a divergent section to the nozzle; a divergent section 
must be fitted instead of a convergent section (which would only reduce the mass 
flow) because when a supersonic stream is accelerated the rate of change of velocity 
is less than the rate of change of specific volume and so to preserve the same mass 
flow the streamlines must diverge. If the ratio of the exit area to the minimum (or 
throat) area is correct for a given pressure ratio, then the gas will issue from the 
nozzle into the receiver in a uniform, supersonic stream. 


The flow of a gas through a convergent-divergent nozzle undergoes more 
changes as the pressure ratio is altered than that through a convergent nozzle. The 
photographs of Figs. 1-4 (pp. 351, 352) show the flow through a convergent- 
divergent nozzle, in which the nozzle has a one-inch diameter throat followed by a 
straight-sided divergent cone with an exit area corresponding to a Mach number of 
1:3. For pressure ratios below the critical the flow is subsonic throughout (Fig. 1) 
and the gas accelerates in the convergent part and decelerates slightly in the 
divergent part. When the pressure ratio is equal to that for which the nozzle was 
designed, as already stated, the gas emerges as a uniform, supersonic jet (Fig. 3). 


If, starting at the design pressure ratio the receiver pressure is increased, then 
ideally, since the flow in the jet is supersonic, the increase in receiver pressure will 
not be felt inside the nozzle and the pressure in the exit plane of the nozzle will 
not change; the nozzle is then said to be in the over-expanded condition. The 
pressure rise which must occur is brought about by a shock wave which springs 
from the nozzle lip, and once again a cyclic flow is set up by reflection. In practice, 
the flow in the nozzle is not everywhere supersonic; at the walls a boundary layer 
exists in which the flow must become subsonic at some level. This allows the 
receiver pressure to be transmitted up the divergent cone and so the shock creeps 
into the nozzle and establishes itself in an equilibrium position along the cone, the 
position depending upon the pressure ratio. The shock is not an oblique one either, 
but frequently consists of two oblique shocks joined by a normal shock. This type 
of flow is illustrated in Fig. 2, where the pressure ratio is such that the shock is very 
near the exit. Notice also that the effects of shock losses and turbulent mixing 
cause the cyclic flow to die out rapidly. 


When the pressure ratio across the nozzle is greater than the design ratio the 
flow in the jet, which is said to be under-expanded, is similar to that from a con- 
vergent nozzle at pressure ratios above the critical (Fig. 4). Again the cyclic flow 
is broken up by mixing. 


In most applications it is not the upstream but the downstream conditions 
which are fixed, the pressure ratio being increased by increasing the upstream 
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stagnation pressure. The descriptions just given apply to this case, except that here 
it is possible to increase the mass flow by further increasing the pressure ratio after 
the critical ratio has been reached, in fact the mass flow then becomes proportional 
to /(P,p,), where P, and p, are the free stream total pressure and density respec- 
tively. The word “choked ” now loses most of its meaning and it is better to talk 
about “ under-expanded,” particularly as this gives a neat sequence of words for the 
convergent-divergent nozzle, namely over-, fully-, and under-expanded. 


1.2. NOISE CHARACTERISTICS 


[Note: In acoustics, sound intensities are normally expressed as ratios to a 
standard intensity which is roughly at the lower threshold of hearing; it is the 
intensity of a sound wave in which the root mean square pressure amplitude is 
0:0002 dynes/cm.”. The response of a human being to stimuli such as sound, light, 
pain, and so on, is roughly logarithmic, and so in measurements it is customary 
to use a logarithmic scale. For the measurement of sound the unit of convenient 
size is the decibel (dB.) which is defined by 


Intensity in dB.=10 log,, (7/I,), where I, is the standard intensity.] 


The noise from a subsonic jet has no discrete frequency components, and when 
it is analysed with respect to frequency there results a spectrum which is fairly flat 
over a considerable range. The slight peak which does occur is at a frequency which 
depends upon the velocity and diameter of the jet; it increases with velocity and 
decreases with diameter and subjectively the jet produces a roar or a hiss, depending 
upon these factors. 


The mechanism which produces the noise was first described by Lighthill, who 
recognised that, in a turbulent flow, that part of the fluctuations of momentum flux 
through a small volume of space which would not occur in a medium at rest and 
transmitting sound constitute, in effect, an externally applied forcing stress system. 
Using this basic idea Lighthill showed that, theoretically, the noise could be 
calculated. Such a calculation involves information on flow patterns which is not 
at present available, but by an approximation he demonstrated certain facts which 
are admirably borne out by experiment, principally that the acoustic power 
depends upon the velocity of the jet to about the eighth power and that the 
maximum sound is radiated at an acute angle to the jet direction. 


But it is not primarily with this sort of noise that this paper is concerned. 

When the pressure ratio across a convergent nozzle rises above the critical 
value, just as the flow pattern changes, so does the noise generated. It is fairly 
certain that the rapid rise in noise level which occurs is due to two effects. 


The first of these is the interaction of turbulence with the shock waves in 
the flow. This phenomenon has received quite a fair amount of attention 
theoretically and is currently being investigated experimentally at Southampton 
University’, where, using a Toepler schlieren set-up on a shock tube, the sound 
waves resulting from the interaction of a single vortex and a shock wave have been 
successfully photographed. 
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NOISE MEASUREMENTS 


The second has a mechanism essentially similar to the edge-tone. Sound waves 


' which arise from the interaction of flow disturbances with shock waves downstream 


are propagated towards the jet orifice and, on arriving there, affect slightly the 
pressure ratio. This can give rise to a small disturbance of the flow and so a feed- 
back loop exists. When the conditions are right, the system can become self- 


' stabilised and the result is that at certain pressure ratios the jet gives out a piercing 


and most intense note of a particular (though not necessarily steady) frequency. This 
mechanism has been investigated by Powell. 


Considerable attention has been given to the reduction of the noise from under- 
expanded convergent nozzles and the results are quite encouraging. They all 
depend to a greater or lesser extent on the facts that the noise generated is sensitive 
to the conditions at the exit of the nozzle (for example a nozzle with a slightly 
crumbled exit lip will be appreciably quieter than one with a sharp, clean edge), 
and also to the rate of mixing, which affects the number of shocks. The devices 
which have been tested result in a jet whose thrust is about the same as, or less than, 
the original one, and up to 12 dB. quieter. 


In this last sentence lies a strong argument for convergent-divergent nozzles, 
for besides eliminating shocks at a particular pressure ratio the nozzle will give 
increased thrust. There is here potentially a considerable advantage, provided that 
the noise from the convergent-divergent nozzle at pressure ratios other than its 
design ratio is, at most, no worse than the convergent nozzle. It is with this problem 
that the present paper is concerned. 


2. Experimental Arrangements 


The general set-up at Southampton University has been described elsewhere 
and present remarks will be confined to the nozzles, measuring gear and running 
technique. 


At present, the air storage capacity limits the pressure ratio available to a value 
corresponding to a (fully-expanded) Mach number M of about 1-5; at this pressure 
ratio the running time for a 2 in. diameter jet is as short as can be conveniently 
handled. Thus in the present tests five nozzles, giving exit Mach numbers of 1:1, 
1:2, 1:3, 1-4, 1-5, have been tested. It was proposed originally to design these 
using the tables provided in Ref. 3, which have been calculated using the method 
developed by Foelsch, but after a little reflection it was clear that with an exit 
diameter of 2 in. it was not possible to work to the required degree of accuracy. 
Accordingly a compromise was effected and the nozzles were made with a short 
conical length of the correct cone angle immediately after the throat, followed by a 
smooth curve into the exit diameter. 


The noise measurements were made with an apparatus made by Briiel and 
Kjaer consisting of a small high-stability capacitor microphone, a third-octave band- 
width spectrometer with 27 filters centred on frequencies from 40 c./s. to 
16,000 c./s., three weighted networks and a linear position, and a level recorder. 
For about half of the work the spectrometer was set to the linear position and used 
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as an additional microphone amplifier to give “overall” response. This involves 
certain inaccuracies (which are discussed in the Appendix) which may be resolved 
by taking spectra. For spectrum measurements the level recorder can be mechanic- 
ally coupled to the scanning mechanism of the spectrometer such that a spectrum 
(to which microphone corrections must be applied) is plotted directly on the 
recording paper. 

For general runs, with the spectrometer to “linear,” a blow-down technique 
was used. The nozzle was taken up to the maximum forward pressure and then 
the pressure gradually reduced. At suitable intervals readings were taken. This 
method is quite useful provided that the pressure is not reduced too quickly, 
otherwise errors creep in owing to reverberation. This point was checked. 


3. Results 

Because of the limited size of the laboratory, a complete examination of the 
noise fields is not practicable, but it is possible to infer a great deal from tests on 
jets where measurements are taken at 30° and 90° to the jet direction. These angles 
are chosen because the maximum noise from a subsonic jet is radiated at (about) 
30° and with an under-expanded jet there may be a strong radiation at 90°. A 
standard radius of 4 ft. has been used; at this distance reflection errors are thought 
to be small. 

Figure 5 shows the characteristics at these positions of the convergent nozzle 
and the convergent-divergent nozzle designed for a forward pressure of 32:5 Ib. /in.° 

The lines for the convergent nozzle display the usual features, namely a rapid 
increase of noise after the critical pressure ratio is passed. Notice that this is 
particularly true of the 90° position; the reason for this is that, with noise associated 
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OVERALL NOISE LEVEL. dB RELATIVE TO 0-0002 


125 
120 
ALL AT 4FT, 30° 
15 
10 20 30 40 50 


UPSTREAM STAGNATION PRESSURE. LB oC) 


Fig. 6. , 
Overall noise level. Several convergent-divergent nozzles designed for different working 
pressu.e ratios. 


with the shock formations (particularly the resonant system which produces the 
whistles), a strong radiation occurs at 90° to the jet direction. The repeatability of 
readings in the under-expanded condition is usually poor when a blow-down 
technique is used because the whistles, which can raise the overall noise level by 3 
or 4 dB., are quite sensitive to pressure and so may be missed in runs of this sort, 
where readings have to be taken at intervals while the pressure is falling. 


The results for the convergent-divergent nozzle are interesting. It will be seen 
that the noise level relative to that for the convergent nozzles gradually decreases to 
a minimum around the design pressure and then rises rapidly. At the design 
pressure the sound at the 30° position is about 6-7 dB. down; at the 90° position 
the improvement is even more spectacular, 10-12 dB. down. (This great improve- 
ment at the 90° position is important with regard to adjacent structures.) The 
comparative performance of the two nozzles is quite well represented here and 


_ apparently there is considerable advantage to be gained with a convergent- 
_ divergent nozzle. 


Figure 6 shows similar results, this time only for the 30° position, for the 
complete set of nozzles. Except for the M=1-1 nozzle, the noise from each falls 


_ 0 a minimum around the correct pressure ratio. It is understandable that the 
| M=1-1 nozzle does not perform quite so well (although it has a minimum, it 


appears at the wrong pressure), because the area ratio is so small and the nozzle 
therefore is difficult to make accurately. 
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Fig. 7. 
Spectrum levels of the noise from a convergent jet. Measurements taken at a point 4 ft. from 
the nozzle exit along a line at 30° to the jet direction. 


Before these results can be accepted it is necessary to make spectrum measure- 
ments to ensure that the different nozzles are comparable on an_ overall 
measurement. Spectrum measurements will also bring out more clearly the 
behaviour of the flow, and so assist in the explanation of the improved behaviour of 
the convergent-divergent nozzle. 


To deal first with the convergent nozzle, five spectra taken at the 4 ft., 30° 
position for different pressure ratios are given in Fig. 7. These were measured with 
the third octave bandwidth analyser and have been corrected for microphone 
characteristics and bandwidth. (The reason for the bandwidth correction is given 
in the Appendix.) It is seen that all the spectra are basically fairly flat, but for the 
under-expanded cases certain discrete tones are present. These latter arise from 
the fact that the resonance mechanism mentioned earlier is working at thes 
pressures, i.e. the jet is emitting a loud whistle. The effect of these discrete tones, 
being as they are at frequencies where the microphone characteristic is beginning 10 
droop, is to make the overall measurements not strictly comparable with those it 
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the range of pressures below about 15 lb./in.?, where the whistles do not occur. 
Because the microphone response is down slightly, full weight is not given to the 
signals in this region. However, this is not to say that the overall graphs are 
meaningless, but that the sound in the upper pressure ranges is probably slightly 
underestimated, relative to that at the lower pressures. 


In Figs. 8 and 9 the convergent and convergent-divergent nozzles are 
compared; in Fig. 8 for the 4 ft., 30° position and in Fig. 9 for the 4 ft., 90° position 
At 38 lb./in.?, there is not much to choose between the two nozzles, but at the 
design pressure of 32 lb./in.? the convergent-divergent nozzle is clearly better. The 
convergent nozzle is here whistling strongly, as can be seen from the two distinct 
tones at the 30° position, and the radiation at 90° is very strong. On the other 
hand, the convergent-divergent nozzle is not whistling and as a result the spectra are 
fairly flat. From the remaining spectra it appears that the convergent-divergent 
nozzle has the advantage for pressures between the design pressure and the critical 
although it has already been seen that the amount varies. Again note that, where 
the spectra for the convergent nozzle have discrete frequency components and those 
for the convergent-divergent nozzle have not, the difference between the overall 
sound levels will be underestimated slightly. 


4, Discussion 


The tests reported here have shown that, for a cold jet operating at a particular 
pressure ratio, there is a definite advantage in using a convergent-divergent nozzle 
designed for that pressure ratio instead of a simple convergent nozzle, certainly for 
exit Mach numbers greater than about 1:3. To take a typical result, a nozzle 
designed for an exit Mach number of 1-4 results, at its design pressure, in an overall 
intensity 6 or 7 dB. less at 30°, and 10 or 11 dB. less at 90° than when a convergent 
nozzle is used at the same pressure. Further, for pressure ratios below the design 
ratio there is still a considerable advantage to be gained, i.e. the provision of a 
convergent-divergent nozzle designed for a particular pressure ratio does not result 
in more, but less, noise as the pressure ratio is built up to the design ratio. 


The extension of these ideas to large-scale hot jets, for example those from jet 
engines, must be viewed with some caution because experimental data are lacking 
for the under-expanded case. In that work which has been done, the whistle 
associated with small, cold jets has not yet been found and since the velocities in 
such jets are very much higher and the “ background ” (i.e. turbulent mixing) noise 
tises very quickly with increase in velocity, noise due to turbulence-shock inter- 
action in such cases may not be of such importance. However, with hot jets there 
may well be another reason for avoiding the formation of shocks, for it has been 
shown by Powell that temperature fluctuations can interact with shock waves to 
produce noise. 


Greatrex® has reported tests on jet engines where the jet was under-expanded, 
although the pressure ratio is not very much greater than the critical. His results 
show that the data in the subsonic range tie up with those in the under-expanded 
case if the latter are plotted at the fully-expanded velocity. The maximum pressure 
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ratio he was able to work at was 2:25, corresponding to a forward pressure of about 
18 lb./in.” in the present work. Figure 5 shows that this is rather low to attempt to 
draw any definite conclusions. Greatrex also tested a convergent-divergent nozzle 
at this pressure ratio and found negligible differences from tests with a standard 
nozzle, but the pressure ratio corresponds to an exit Mach number of only 1:14 and 
it is difficult to get reasonable results at these Mach numbers. 


The question of the applicability of the present work to full-scale jets is an 
important one, for jet engines are likely to work at pressure ratios of the order 5:0 
(exit Mach number 1°7) in the not too distant future, and the reductions offered are 
certainly worth while. The point could easily be settled by some tests on model 
hot jets. 
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Appendix 


The overall response of a microphone and amplifier depends upon the frequency 
characteristics of both the microphone and the amplifier and upon the shape and 
frequency range of the noise spectrum under observation. Thus, in Fig. 10(a), the 


Fig. 10, 
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overall response measured with a microphone and amplifier whose joint frequency 
response is as shown by curve A will be different for the two spectra B and C, even 
though they may contain the same energy. For this reason, overall measurements should 
be treated as comparative and not as absolute ones, particularly where discrete tones are 
present, e.g. in Fig. 10(b) the discrete tone in spectrum D will not be weighted correctly 
and so the difference between the two spectra will be underestimated. 


In the present work, because the spectrum does not move much, in the subsonic 
range the readings are comparable to less than 1dB. The readings in the under-expanded 
range may be underestimated by as much as 2dB. relative to the subsonic levels. 


This difficulty does not arise with spectrum measurements, for then the (microphone 
and filter) corrections may be applied and results given whose absolute accuracy depends 
only upon the absolute accuracy of the overall calibration; but the result measured 
in a band depends upon the width of band and if results from different analysers are 
to be compared, corrections should be made for bandwidth. This has been done 
in the present work, levels having been reduced to spectrum levels using the usual 
expression 


Measured level in filter of 
bandwidth Af cycles 


To compare results given here with those from any other analyser, the appropriate 
value of 10 log,, (Af) should be added. 


Spectrum level = [ | — 10 log,, (Af). 
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Moments and Deflections of a 
Simply-Supported Beam Grillage 


J. P. ELLINGTON and H. McCALLION 


(Departments of Civil and Mechanical Engineering, University of Nottingham) 


SUMMARY: By using the methods of the Calculus of Finite Differences, 

expressions are obtained for the nodal moments and deflections of a simply- 

supported grillage, subjected to a loading constant along one set of beams and 

having a sinusoidal variation along the other set of beams. A simple example 
verifies the expressions and illustrates their use. 


1. Introduction 


The use of latticed beams or grillages in engineering practice is well known; 
they are found in bridges, foundations, structures, ships’ hulls and aircraft fuselages. 
The problem of calculating the moments and deflections of the component beams is 
difficult, however, and generally it is necessary to use either numerical methods" 
or approximate methods. 


Hetenyi® has recorded investigations by himself and Pippard, while recently 
Morice has discussed and extended work by Guyon and Massonnet. In general 
these approximate methods take one of three forms. One set of beams may be 
replaced by a uniform continuous medium whose properties simulate the original 
discrete beams, or the loading due to a set of beams is replaced by a series of 
continuous sinusoidal loads. Hendry and Jaeger, have combined these methods 
and obtained a convenient design procedure by assuming that the deflection and 
loading on the longitudinal beams follows a simple sinusoidal curve. Timoshenko’ 
has also established that both sets of beams may be replaced by a continuous 
medium, and the the gridwork may be considered as an aleotropic plate. 


However, to the authors’ knowledge, no exact solution has been obtained for 
the deflections and moments of a grillage, and the purpose of this paper is to present 
such a solution for a particular case. The grillage examined is taken as being 
made up of two orthogonal sets of equally spaced uniform beams, attached at the 
nodes in such a manner that, while the beams never part, they do not transmit 
bending or torsional moments to each other. The loading is taken as a set of 
concentrated loads applied at the nodes, the loading being constant along one set 
of beams but having a sinusoidal variation along the other set, summation of sine 
terms allowing a wide range of loading conditions. 


Received March 1957 


360 The Aeronautical Quarterly 


as u 

then 

NoT. 

2. 

be 

in 
i 

of 


n) 


erly 


SIMPLY-SUPPORTED BEAM GRILLAGES 


The moments and deflections at the nodes of the two sets of beams are taken 
as unknowns, and the two equations of continuity and an equilibrium equation are 
then solved explicitly, using the methods of the Calculus of Finite Differences. 


NOTATION 
(r,s) nodal point of orthogonal sets of beams in r- and s-directions 
R,S total number of r-beams is (R - 1) and of s-beams is (25-1) 
I,,I, second moments of area of r- and s-beams 
L,,L, spacing of beams in s- and r-directions 
E Young’s modulus 
nodal deflection 
M,,, nodal bending moment in r-beams 
m,,, nodal bending moment in s-beams 
P_ concentrated load at the central nodal point 
W,,., nodal concentrated load 
w assumed form of loading, W,,=w sin (mtr/R) 
n an integer 


E,,E, finite difference “ shifting ” operators, defined’ by 
E,*E.¢y, r,s Vr+a,s+z 


4. differencing operator defined by A=(E-1) 
¢ finite difference operator, (E- 1)?/E 

k flexural stiffness, 12E7/L* 

A= (k,/k,) {[cos (nz/R)- 1]?/[2 cos (nz/R)+4]} 


she B| disposable constants in general solutions of difference equations 


characteristic exponents 
a,8 characteristic exponents, 7, ,=2+i68 
a,b constants, a=(1-A), b= /(6A—A?) 
y=arce tan (b/A) 
characteristic constant, =(6E/,/L,”) [(cosh 1)/(cosh o + 2)] 


2. The Difference Equations 


Consider now a rectangular grillage of simply-supported uniform beams. The 
beams in the r-direction are of modulus of elasticity E and second moment of area 


_ 1, and are spaced at equal distances apart L,. These beams rest on cross-beams 


in the s-direction, at right angles to them, the cross-beams having a second moment 
of area J,, and being spaced L, apart. 
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The deflection at a point where the beams cross, (r, s), is y,,,, it being assumed 
that the two sets of beams are simply attached to each other without any torsional 
or bending restraints. The bending moment in the r-beams at the intersection (r, s) 
is M,,, and the moment in the s-beams is m,,. Then, if the loading is taken as 
being applied only at the nodes of the grillage, Clapeyron’s three-moments 
equation“ may be written for two adjacent spans of a beam in the r-direction as 


M,+2,s M,,.= (6EI,/L,”) 2Yr+1,8 2). 
Similarly, in the s-direction, 
My, s+2 4m, 541 (6EI,/L,”) y + yr, 


An equilibrium equation can then be written for the shear forces and applied 
loading at a node (r+1, s+1), 


041 2M 41,041 M,, Mr+1, — 241,041 + Mr+1,8 


W,., being the concentrated load applied at (r,s). 


These equations are a set of non-homogeneous linear partial difference 
equations” in the unknown deflections and bending moments and, using the partial 
finite difference “shifting” operators E, and E, defined by E,*y,.=y,42,, and 
E,*y,,s=Yr,s+2, they may be written in operational form as 


L,’ 2 L, 2 one 
1)?m,,.+ E —(E, - 1}*y,,.=0 (2) 
and 1)? E,M,.,.+ . (3) 


Since there is now no danger of ambiguity, the suffixes r, s will be dropped 
from the symbols for deflections, moments and loads. Equations (1) and (2) can 
then be written as 


_ 6El, 


Substituting equations (4) and (5) in equation (3) gives 
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12EI (E- 1)? 
Here k= and 
It is now necessary to know the form of the applied loading, and it will be 
assumed that the loading along an s-beam is constant and has a sinusoidal variation 
along the r-beams. Obviously a wide range of loading conditions can be repre- 
sented by the sum of a series of sine terms, but for convenience only one term is 


used below and it is assumed that 
W =w sin (ntr/R), ‘ (8) 


Here (R - 1) is the number of cross-beams and n is an integer. 


(7) 


The partial difference equation (6) may now be solved by taking 


where y, is a particuar solution satisfying the non-homogeneous equation and the 
boundary conditions in the r-direction, while y, is a complementary solution satis- 
fying the homogeneous equation, the sum of these solutions being subject to the 
boundary conditions of simple supports, or zero deflections and moments at the 
edges. 


Suppose now that y,=Y sin(mtr/R)= exp(intr/R), . (10) 


where Y is to be determined from equation (6) 


{= (5) + Y exp (inzr/R)=w exp (inzr/R). 


Using now the property of the E operator that f (E,) a’ =a’f (a), 


exp (intr/R)= (2 cos (nz/R)+4] 


exp (intr/R). 
To solve the second operational term, E, is replaced by (1+A,), 4, being the 
partial difference operator defined by — 


exp (inzr/ R)= (A, + 1) (4,2 + 64, +6) 


Then exp (intr/R). 


If this operational expression is then expanded in ascending powers of A,, the 
first term is A,*/6 and, since it is a partial operator, it, and succeeding terms, make 
no contribution. 


The value of Y is then found from equation (11) and 


W [cos (nz/R)+2] 
[cos (nz/R)— 1}? 


sin (ntr/R). ‘ (12) 


This obviously satisfies the conditions y,=0 at r=0 and R, and also 
M,=¢ (E,) y, =0 at r=0 and R. 
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To solve the homogeneous equation 


{ & +e } 


assume that y,=A exp (cs) sin (ntr/R), . (14) 


where o is to be determined, and A is a disposable constant. Equation (13) then 
becomes 


"(2 cos (nz/R)+4] *(2 cosh +4) 


[cos (nx/R)— 1}? 


Writing k, [2 cos (nz/R)+4] 


cosh’o —2 (1 - A) cosh o +(1+4A)=0 
or cosh -A)+'/(A?— 6A). ‘ (16) 


This is a characteristic equation giving a total of four roots, there being a 
positive and negative value of « associated with each coshc. 


It can be seen that three forms of solution are possible, depending on the value 
of A. When A <6, as will usually be the case unless k,/k, is disproportionately 
large, the roots of equation (16) are complex; when A > 6 the roots are real and 
unequal, and when A=6 the roots are equal. 


3. Complex Roots, A < 6 


Taking the most practical case of A < 6, giving complex roots, the character- 
istic equation (16) becomes 


cosh ,=a@+ib, 
where a=(1-A), b= / (6A—A?) and i=/J(-1). . 


The corresponding values of o are then +(2+i8) and +(a-if), where 
cosh «cos B=a and sinh z sin B=b. 
Squaring and adding these equations, 
cosh?a + cos?B= a? + b? +1. 
Then, adding and subtracting alternately 2 cosh «cos 8=2a gives 


2 cosh a= {(a+ 1)? +b?} + {(a- 1)? +57} 
= /(44+2A)+ / (6A) 
and 2cos B= /(4+2A)- / (6A) 


364 The Aeronautical Quarterly 


(18) 


an 


fol 


th 


CO 


eq 
T 
ir 
T 

a 


5) 


y 
d 


SIMPLY-SUPPORTED BEAM GRILLAGES 


Returning now to equation (14), 


y.=sin (ntr/R) > A exp (cs) 
=sin (ntr/R) {exp (as) [A cos 8s +B sin Bs}+ exp (- 2s) [C cos Bs+D sin Bs]} (19) 


The total deflection y is then the sum of equations (12) and (19). From 
equation (4) the bending moment in the r-direction is 


_ 6EI, [cos (nz/R)- 1) _ [cos (nt/R)- 1] 
L,? [cos (nz/R)+2] [cos (nz/R)+2} 


and, from the form of y,,,, the boundary conditions in the r-direction are satisfied. 


M (20) 


An expression for the bending moment in the s-direction may be developed as 
follows, from equation (5), 


E, - 1)? 
{ } (y, + y.) =sin (ntr/R) (cs), 


_6El, (cosh 1 
\cosho+2/’ 


there being two values of » corresponding to the two values of cosho. Using 
equation (17) gives 


where 


= /6) ib)=(k.L,/6) (6X) exp 
where y=arc tan (b/A). 


Thus m,,,=k,L,/ (A/6) sin (ntr/R) {exp (as) [A cos (8s+y)+B sin (8s+y)] + 
+ exp (- 2s) [C cos (8s—y)+D sin (Bs—y)]}. . (21) 

It is now convenient to consider a grillage with an odd number of cross-beams 

in the s-direction; the origin of s is then taken at the centre, the supports being at 
s=+S. Symmetry then demands that the deflections and moments on either side 


of the origin should be equal or y,,,=y,,-, and m,,,=m,-,, giving A=C and 
B= -D. 


Thus y,,,=2 sin (ntr/R) { +A cosh zs cos Bs +B sinh as sin Bs } 


The zero moment condition at the support gives 


A ___tanh aS tan + tany 


1—tany tanh 2S tan BS’ 
and the zero deflection condition gives 


w sech a§ sec BS cot y (tan y tanh 2S§ tan BS— 1) 
4k,A (1 + tanh?2S tan?) 


__wsech 2§ sec BS cot y (tanh aS tan BS + tan y) 
4k,A (1 + tanh?2S tan?8S) 
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4. Real Roots, A > 6 


When the value of A is greater than six, the characteristic equation shows that 
both values of cosh o are negative and it is necessary to use the transformation 


exp (c)= — exp (8). 
The characteristic equation is then 
cosh 6, ,=(A—1)F /V(A?-6A). , (23) 


Again, taking the origin for s at the centre, the symmetry condition allows only 
a solution containing even functions, so that 


sin (ntr/R) +(-1)[A cosh 6,5+ B cosh 6,5] } 


and m,,.= sin (ntr/R)(- {u, A cosh 6,5+4,B cosh & 


6EI, ( 
where b= - 


1+cosh ;) 
2-cosh 6/° 


The corresponding values of A and B are found from the support conditions as 


cosh 6,5 4k,A 


In the unlikely case of A=6, the values of @ obtained in equation (23) are 
equal, thus yielding only two solutions instead of the required four. Consequently, 
as in the theory of linear differential equations, the original assumed solution, 
equation (14) must be replaced by 


y,=(— 1} (A + Bs) exp (@s) sin (ntr/R), : ; . (26) 
and a solution can be developed as before. 


5. Application of Method 


The use of the equations, and also their validity, is best demonstrated by a 
very simple example whose solution may be checked by ordinary beam theory. 
Consider two orthogonal beams, simply-supported at their ends, and carrying a 
concentrated load P at the central nodal point. The lengths of these beams will be 
2L, and 2L, respectively and the flexural stiffness, (k=12EI/L*), of the beam in 
the r-direction will be taken as twice the stiffness of the beam in the s-direction. 


Thus, using the notation of the paper, R=2, S=1, the load P being at 
(r, s)=(1, 0), the supports of the r-beam being at (r, s)=(0,0) and (2,0), while the 
supports of the s-beam are at (r, s)=(1, 1) and (1,- 1). 
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SIMPLY-SUPPORTED BEAM GRILLAGES 


From equation (8), the sinusoidal variation of the load is given by 
W,..=w sin (ntr/R)=P sin (r/2). 


In more complex problems it will be necessary to take several sine terms to 
represent the actual loading. 


Remembering that k,=2k,, the parameter A is now calculated from 
equation (12):— 


__k, [cos (nz/R)- 17 
k, [2 cos (nt/R)+ 4} 


The characteristic exponents 2 and @ are then found from equation (18) as 


VOL 
The angle y is then taken from equations (17) and (21) as 


tan y=b/A= [(6/A)- (1D). 


The two disposable constants A and B in the expressions for moments and 
deflections are then found, following equation (22), as 


w sech aS sec BS cot y (tanh 2S tan BS'+ tan y) 
4k, (1+ tanh?aS 


w sech 2§ sec BS cot y (tanh 2S tan BS tan y - 1) 
4k,A (1 + tanh? tan?@S) 


[SP/ {6k,/ (11)}]. 


From equation (22), the nodal deflection at a junction (r,s) is 


yr,s= 2 sin (ntr/R) {(w/4k,A)+ A cosh as cos Bs +B sinh as sin Bs}, 


and so in this case the deflection at the junction of the two beams (r, s)=(1, 0) is 


P 2P 
Vi,0=2 +A] 3k, 


The moments in the r-beam are found from equation (20) as 


_ L,k, [cos (nz/R)-—1] 
[eos (nz /R) +2] 


Thus the bending moment at the centre is 
—(L-k,/4) y,,5=—(PL,/3). 


23 
nly 
24) 
25) 
= — (P/6k,), 
= - 
ry. on PL,’ 
= 
in 
at 
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The expression (21) for the bending moment in the s-beam may be simplified to 
m,,,=(Ak,L,/3) sin (ntr/R) {A (cosh as cos Bs - tan y sinh as sin Bs) + 
+B (sinh as sin Bs + tan y cosh as cos Bs)}, 
and substituting (r, s)=(1, 0) gives 

m,,,=(Ak,L,/3)(A +B tan y) = —(PL,/6). 
These solutions are readily checked by straightforward beam theory, and 
indicate the ease with which an exact solution may be found. It is realised that 
the method is not a general design method, since not all possible loading conditions 


are considered, but the method is capable of extension and could be of great use 
in considering the validity and accuracy of existing approximations to grillages. 
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On the Generation of the Streamwise 
Component of Vorticity for Flows in 
Rotating Passages 


AUSTIN G. SMITH 
(Imperial College, London) 


SUMMARY: A method is presented for the calculation of the streamwise 
component of vorticity, for flows in rotating passages. The method may be 
regarded as an extension of the methods applied in recent years to the 
calculation of secondary flows in stationary passages. 


Attention is concentrated on the quantity ]=p+4pV? — 4pU?, where V is 
the fluid velocity relative to the rotor and U is the rotor tangential velocity. 
I remains constant along a streamline in the rotor, and is found to enter into 
the equation for the generation of vorticity in much the same way as the total 
head enters for flow in stationary passages. 


An approximate calculation of the streamwise vorticity generated in a 
simple axial-flow rotor is made, and qualitative consideration is given to the 
flow in a centrifugal impeller. Just as for the calculation of secondary flows 
in stationary passages, an approximate shape of the streamline must be 

assumed before the secondary flows can be calculated. 


1. Introduction 


Experimental studies of the flow of fluid in curved passages usually show that 
near the walls the fluid direction deviates markedly from that in the “ main stream.” 
These deviations are often explicable qualitatively in terms of the smaller transverse 
pressure gradient necessary to turn lower-velocity elements of the flow, or in terms 
of the “ centrifugal force ” on low-velocity elements. The first important quantita- 
tive study was that of Squire and Winter who showed that by consideration of the 
component of vorticity in the flow direction (termed hereafter the “ streamwise 
component ” of vorticity) it was possible to calculate secondary flows in a stationary 


cascade of “impulse” blades. The secondary velocities were calculated from the 
streamwise component of vorticity by defining a stream function ¥ for the secondary 
Velocities and solving V?¥=—wy where wy is the streamwise component of 


vorticity, for the passage formed by the blades. Experimental results confirmed 


their theory. 


Hawthorne™ produced a more general method for the calculation of streamwise 
vorticity and, deriving largely from his work, a better understanding of the nature of 


4 Received July 1956 


November 1957 369 


173, 


AUSTIN G. SMITH 


secondary flows in passages and in cascades of aerofoils has been achieved“: “. The 
cascade theory uses ideas from actuator disc theory and a consideration of the 
strength of the vortex sheet shed from the blade trailing edges, as well as the stream- 
wise vorticity theory. The success of the theories is not complete, since full 
explanation of the losses of total pressure which appear to be associated with 
secondary flows is not achieved and, further, the theories essentially ignore viscous 
and turbulent stresses. However, it appears that some advance in the understanding 
of flows in rotating passages will be possible by extension of the methods which are 
at present applied to stationary passages. 


The first step is to find a method for the calculation of the streamwise vorticity 
for flow through rotors. The present paper gives a method for such calculation. 


NOTATION 


a angle between the relative streamline, or an arbitrary line, and the 
generator of the cylindrical surface containing the line (Fig. 3) 


C length of a circuit 


y angle between the direction of curvature of the relative streamline 
and (V//p), V (Fig. 1(@)) 


8 angle between @ and VI/p 
6 polar co-ordinate (Fig. 3) 


D a 
Di operator a t V.V or at v.V 
z partial differentiation with respect to time 


H total pressure (=p + 4pv?) 
I quantity defined as p+4pV? -4pU? 
i,j,k unit vectors in the direction of the x, y and z axes respectively 
a angular velocity of the rotating co-ordinates; scalar, 2 
static pressure 
density of fluid 


position vector 

radial distance from the z-axis (Fig. 3) 
principal radius of curvature 

distance along a streamline or arbitrary line 
circumferential velocity (=@2,r); scalar, U 
velocity in the rotating co-ordinates; scalar V 
velocity in the stationary co-ordinates; scalar, v 


Za dau Be 


vorticity in the rotating co-ordinates; scalar, W 
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w vorticity in the stationary co-ordinates; scalar, w 


x,y,z Cartesian co-ordinates (Fig. 3). 


| Suffixes 


1,2 at position 1 or 2 
6,r,z in direction of increasing 4, r, z 
V in direction of the relative velocity V 
vin direction of the absolute velocity v 
secondary 
‘,y,Z in direction x, y, or Z. 


2. The Equation of Motion in Terms of the Quantity I 


An essential feature of Hawthorne’s theory” is that the gradient of the total 
pressure is related to the growth of the streamwise component of vorticity. The 
fact that total pressure is constant along a streamline for inviscid flow makes it 
possiblet to treat this relationship. In rotating passages, such as those between 
turbine or compressor blades, the total pressure is not constant along a streamline. 
In order to find a relation similar to Hawthorne’s it is necessary first to put the 


- equation of motion into a form containing a quantity which is constant along a 


streamline passing through the rotor. This is done in what follows, the resulting 
equation being equation (9). 


For inviscid, incompressible flow in the absence of body forces the equation of 
motion is 


Dv 1 


Taking a co-ordinate system rotating with steady angular velocity 42 and 
instantaneously coincident with the stationary system, then the acceleration Dv/Dr 
of a particle in stationary co-ordinates is related to its acceleration DV/Dt in the 
rotating co-ordinates by 


Dv DV 
Di Dr +22,.V+Qa (<2 (2) 
The equation of motion in the rotating co-ordinates is therefore, combining (1) 
and (2), 


1 DV 


| There are various relations by which (3) may be thrown into different forms. 
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The relations are:— 


D_a 
The definition Dr= +(V.V) 
The relations for rotating co-ordinates 
2 
The vector identity (V.V) V=V -VaulV. . . . @ 
and the definition w= curl v. 


Manipulating (3) with the aid of these relations, we may obtain 
1 ov 


Equation (9) is the fundamental form of the equation of motion which we 
shall use to develop the secondary-flow theory. The similarity between (9) and 


which is the equation of motion for stationary co-ordinates, is at once obvious. 
Equation (9) reduces to (10) when 42 =0. 


If we now consider flow through a row of blades rotating with angular velocity 
a2, and we suppose that there are no stationary blades near enough to the rotor to 
impose cyclic changes in V, then dV/d0t=0 and the equation of motion is 


Thus the vector product of the velocity in the rotating co-ordinates, and the 
vorticity in the stationary co-ordinates, is equal to the gradient of I/p. The 
similarity between (11) and 


. 


leads us to expect that an equation for the growth of the streamwise component of 
vorticity in rotating co-ordinates may be found by much the same method as 
Hawthorne” used for the case of stationary co-ordinates. 


The quantity 7 is constant along a streamline for steady flow in rotating 
co-ordinates. This follows from (11), since for steady flow D/Dt=V.V and so 


(4) =v.(v-) =V.(Vaw)=0. 
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D/(H\_D 
and hence Dr 
4) By integration along a streamline we obtain the familiar equation (for 
steady flow) 
5) H, - H;=p [((Uvs), (Uv),] ‘ 
6) for the change in total pressure between two points on a streamline in the 


rotating co-ordinates. 


3. Vector Equation for the Growth of Streamwise Vorticity in 
Rotating Co-ordinates 


Equation (11) implies, by the definition of a vector product, that vortex lines in 
3) the stationary co-ordinates and streamlines in the rotating co-ordinates are con- 
tained in surfaces of constant J. The vorticity in stationary co-ordinates (hereafter 
called for shortness the absolute vorticity) may therefore be resolved into 
two components : — 


Je (a) along the velocity in the rotating co-ordinates (hereafter called for 
d shortness the relative velocity), 


(b) normal to the relative velocity but tangential to the surface of constant /. 
There is no component normal to the surface of constant /. 


8. _(V.w)V 


The component (a) is (15) 


V V? 
ly (Va w)V 
Component (5) is 
Now w being curl v, and since div curl=0, it follows that div w=0 and the diver- 
) gence of the sum of the components of w is zero, or 
. By expansion of (15) and a certain amount of manipulation, which is given in 
detail in the Appendix, we may show that 
Wy 2 2 
f 
: This equation, all the terms of which have a simple physical interpretation, is a 
usable vector equation for the rate of change, along a relative streamline, of the 
streamwise component of vorticity. When 42=0, equation (17) reduces to 


This is Hawthorne’s fundamental equation for the growth of the streamwise 
component of vorticity in stationary passages. 
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Fig. 1. 


4. Scalar Equation for the Growth of Vorticity 


Equation (17) may be rewritten as 


Taking the first term on the right hand side of (19), (V.V) V is the fluid 
acceleration at a point on the relative streamline, and it may be resolved into two 
components, one along the streamline and the other in the direction of the principal 
normal n to the streamline. The former may be ignored, as its scalar product with 
(VI/p), V is zero. The magnitude of the component in the direction of the principal 
normal n is V?/R, where R is the principal radius of curvature of the relative stream- 
line. The first term in the right hand side of (19) is therefore 


2 

cosy, 

where y is the angle between the principal normal to the streamline and the direction 
of (VI/p),V, as shown in Fig. l(a). The second term on the right hand side of 
(19) is 


20 | 


where 6 is the angle between V//p and @. Equation (19) may be rewritten 
therefore as 


w 20 


This equation may be aieaial along a relative streamline to give, between points 
(1) and (2), 
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Equation (21) may be used to calculate the course of wy along a relative streamline 
if the shape of the relative streamline and the magnitude and direction V//p 
are known. 


Far upstream of a rotor row, the flow will normally be in “ radial equilibrium,” 
and the surfaces of constant J will be right cylinders. V/J/p will therefore be radial 
in direction, and the angle 6 will be </2, and the second term on the right hand side 
of equation (21) will be zero. This term will remain zero until, under the action of 
the blade row, the surfaces of constant / are distorted from cylindrical. The first term 
on the right hand side of (21) will differ from zero, however, as soon as the relative 
streamline becomes curved by the action of the blade row, independently of any 
distortion of the surface of constant / from the original cylindrical form. 


For approximate calculation of flow through rotor blade rows it may be 
sufficient to assume that the surfaces of constant J are little distorted from their 
original cylindrical shape and, with this assumption, the cosy and R terms can be 
treated. It is interesting to note that if V//p is normal to a rotor blade surface, as 
it will tend to be in the boundary layers on the blades, the value of 5 will differ from 
=/2 in general. When the stagger of the blade is zero, the value of 6 will be near 
t/2. When 6 differs from «/2, the second term in (21) will make a contribution to 
the growth of w;y/V. This is the contribution leading *. the “ centrifuging ” of the 
boundary layer. 


For the arrangement shown in Fig. 1(b), where V.’/p is supposed to be normal 
to the blade surface, 5 is essentially less than 7/2, cos is positive, and positive 
streamwise vorticity is generated by the action of the second term in (21). will be 
z/2, and the first term of (21) therefore makes no contribution to the growth of 
the vorticity. The higher the stagger ¢ of the blade (defined as in Fig. 1(b)) the 
smaller is § and hence the larger the effect of the second term in (21) for a given 
2, V, and VI/p. 


5. Approximate Calculation of the Growth of Streamwise Vorticity 
in a Simple Rotor 


5.1. GENERAL 


Consider a rotor as shown in Fig. 2. The absolute flow is axial at entry to the 
rotor, and the relative flow is axial at exit from the rotor. At entry to the rotor, a 
boundary layer has built up on the feed duct walls such that at the walls the axial 
velocity v, is 0-75 x the main stream axial velocity vn. Away from the walls the 
velocity v, rises linearly to vz at r/r,=0-6 and 0-9. 


The problem is the determination of (w,/V), from equation (21). Considering 
VI/p first, as this is the quantity on which (wy/V), mainly rests, it will be shown 
that VI/p does not change through the rotor, and its value may therefore be 
determined from the approach flow. 


The fundamental method of doing this is use of the equation of motion (11), 
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Fig. 2. 


In the approach flow the only component of w is we, and thus U, w is zero. 
v,W has only the radial component —v,ws and thus 


dv, 
VI/p= 5, (23) 
This result is equally easily derived in this instance from equation (13). Since 


v and U are orthogonal, v. U=0, and hence //p=H/p. 


dv, 
Thus ap 5, 
Regarding the constancy of V//» in passing through the rotor, it is easy to show 
that, for the actuator disc into which the rotor row shrinks when it has an infinite 
number of blades of zero axial length, the velocity distribution v,=0, v,,=v,, at aly 
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radius, and vp,=Or is a solution of the equations of motion and continuity con- 
sistent with the boundary conditions stated in Fig. 2. From v,=0, it follows that a 
streamline does not alter in radial position as between upstream and downstream of 
the disc, and that V//p is the same on either side at a given radius. 


We shall assume this velocity distribution for the blade row with finite number 
and axial length of blades. Implicit in this assumption is that V//p remains radial, 
although in fact the secondary flows in the rotor passages twist the surfaces of 
constant J out of their original cylindrical shape. The validity of this assumption is 
commented on in Section 5.4. 


With V//p radial in direction, the second term on the right hand side of 
2 

equation (21) is zero. We need only evaluate (w,/V), and (FR VI/p | cos y ds 
1 


to determine (w,/V).. 


5.2. VALUE OF (wy/V), 


Wy W.V_ wV,+wVotw, 
We have vy = (24) 
With the co-ordinate system of Fig. 3, the r, 6, z components of vorticity are 

_ 

700 
Ov, Ov, 


Having axial symmetry, 0/00=0. and having radial equilibrium, 0/¢z=0. Thus 
— Ov,/ dr, and w,=w,=0. 


Also = —sin 2,, V./V =cos z,, and (24) becomes 
wr) Ov,, SiN COS 
( 7) Or Vey - (6) 
where a, is defined as in Fig. 2. (w,/V), is thus easily determined, in terms of the 
angle z, commonly used in compressor theory. 


2 


5.3. VALUE OF 


2 
|V 
PR | I/p | cos y ds 
Since it is assumed that the relative streamlines are contained in right 
cylindrical surfaces, the curvature of an arbitrary line drawn on a right cylinder 
must be considered, in order to treat y and R. 
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Fig. 3. 


If ds is an element of arc along the arbitrary line, then with « taken positive in 
the sense shown in Fig. 3, 


dz= ds cos 2 
dx= ds sin sin : ; 127) 
dy= —-dssin zcos 6 

and dr= idx+jdy+k dz. 


Therefore a = isin a sin 6-jsin zcos0+kcosa 


Differentiating (28) with respect to s 


ar ( dé 
sin 2 cos 6 + COs asin 0%) 


4 i( sin a sin 6 d6 COS 
ds ds 


+k (-sin a =). : 


This is the curvature of the line. It may be resolved into two components, one 
radial, and the other normal to the line and tangential to the cylinder on which the 
line is drawn. The radial component, which is 


i( sin a cos 6 =) +5(sin a sin 6 


makes no contribution to the desired integral, since y is +7/2 for this component. 
The other component, for which y is 0 or 7, is 


i(cos « sin 932) cos cos 6%) (sin 
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Its magnitude is da/ds. It might be called the “developed curvature,” since 
it is the curvature which the line would have if the cylinder were developed onto 
a plane. 


The integral is therefore 


32 
+ 2 Vi /p| 


1 


2 
2 
ve 
1 
since v,=V cos a. 


Near the inside diameter of the rotor row, (V//p), V is in opposition to the 
developed curvature, and the growth of wy is negative. The opposite is true near 
the outside diameter of the rotor row. With the assumption that the relative 
streamlines are contained in right cylindrical surfaces, v, is constant, and so is 
|VI/p|. The integral is therefore 


+2 da= | . . (30) 


On the right hand side of (30) the sign has been taken in accordance with 
equation (17). 
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TABLE I 

rir, 0-5 0-55 0-6 0-9 0-95 1-0 

2-5 2:5 2-5 ~2:5 —2:5 

33-6 321 30-9 41-9 47-3 52:9 
~1-35 1-82 2-07 2:37 
Q A(r/rt) 

w, | -0-18 +007 +0°31 


5.4. VALUE OF (wy/V), 


From equations (30), (26) and (21) we have 


2 


(7) = 2 | costa sin dv,/Or 61) 
1 


2 


Uz 


Since 2,=0, v,,=v.. If the rotor speed be set relative to the fluid velocity such 
that Or, =v.m, then (31) becomes 


Wy (ve/ Vem) 
tens 


(32) 


It is W, the relative vorticity, which is of most use for determining the 
secondary velocities. 


Since W=w-2Q, we have Wy=wy —20,=wy - 20 cos 15633) 


and equation (32) becomes, remembering that cos 2,=1, 


(We) 


Values of W,/© calculated from equation (34) are given in Table I. The table 
is carried only over the r/r, values where dv,/0r differs from zero. 


Figure 4 shows graphically the distribution of W,/Q. It will be seen from 
Table I that the wy term is of much the same size as 22. It should be borne in 
mind however that a comparatively small value of 0 (v./v.m)/0 (r/r,) has been taken. 
and that doubling this would make the wy, term predominate. 


Neither the axial length of the rotor nor the circumferential pitching of the 
blades enters the expression for W,. These quantities have an importance, however, 
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in that if the secondary velocities (which for a given distribution of Wy depend 
partly on the circumferential spacing of the blades) are large enough to distort the 
surfaces of constant J substantially from the assumed cylindrical shape, then the 
related important assumption is invalidated. 


An approximate estimate of the rotation of the surfaces of constant J may be 
made by considering the circulation round the circuit ABCD (Fig. 2) in the exit 
plane of the rotor. By Stokes’s theorem 


f V.av= |W. aA, 


A 


where dv is an element of length along the circuit, and dA is an element of the 
surface ABCD. The mean secondary velocity V, around the circuit is therefore 


= W, dA 
(35) 
With the values of W,/Q given in Table I, V, is found to be —0-1lvm. Assuming 
that this velocity applies to the mid-points of the sides AD, BC, then if the blade 
axial length were 0-25r,, the radial movement of a particle reaching the mid-point 
of AD or BC would be about 0-02r, if the radial velocity of the particle increased 
linearly in its passage through the rotor. Movements of this-order would not alter 
VI/p sufficiently to make the calculations leading to Table I substantially in error. 


PROFILE OF AXIAL 
VELOCITY AT ENTRY 
TO IMPELLER 
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6. Qualitative Consideration of the Secondary Flows in the 
Impeller of a Centrifugal Compressor 


The simple rotor considered in the previous section may be taken as the 
“rotating intake guide vanes” or “inducer vanes” of a centrifugal impeller. If 
we imagine that the friction on the feed duct walls has produced a velocity profile 
as shown in Fig. 5, then at point A, near the outer diameter of the impeller channel 
in the intake guide vanes, positive wy is induced by the first term of equation (21), 
while for point B near the inner diameter of the channel, negative wy is induced. In 
this intake region, V//p is approximately at right angles to €2, and the second term 
of (21) is small. 


At point C, on the same relative streamline as A, but in the radial-flow part of 
the channel, then provided the relative streamlines are following paths uniformly 
interpolated between the impeller surfaces, VJ/p is approximately in the same 
direction as €2. The second term in equation (21) will thus contribute positive wy. 
At D the opposite is true. At both C and D, the radius of curvature being large, 
the contribution of the first term in (21) is small. All along the path AC in the 
impeller, positive wy will be generated, while along the path BD, negative wy will be 
generated. The secondary velocities corresponding with these vorticities will tend 
to sweep the fluid with low / from both walls of the inlet duct towards the “leading” 
or “suction” side of the channel. The existence of a region of low J near the suction 
side of the channel has been found experimentally. 


If the process of sweeping the fluid with low J onto the suction side of the 
channel were completed, then V//p would be in the direction of - U and therefore 
normal to €2, so no further contribution of streamwise vorticity from either term in 
(21) would occur, the contribution of the first term being zero on account of the 
zero curvature of the streamline, and that of the second term being zero since 
a .(VI/p)=0. As far as the present theory is concerned, the condition that / 
becomes smaller as the suction side of the channel is approached, and the surfaces 
of constant J be radial, is a stable one. The processes outlined will also result in a 
tendency for the frictional boundary layer generated on the channel walls to be 
swept to the suction side of the channel. 
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Appendix 
DERIVATION OF EQuaTION (17) 


Expanding the terms of equation (15), 


Wy w | 1 
div V+(V.V) + div (Vaw)aV + [(Vaw)V].V =0 (36) 
Equation (36) may be manipulated by means of the relations (37), (38), and (39):— 
divV=0 (continuity) . (37) 


div (V,w),V = V.curl (V,w) — (V,w). curl V 
=V .curl V (I/p) (V,w) . curl (v U) 
Since curl V =0, curly=w and hence (V,w). curl y=0 and curl U=2@2 


1 
V2 


Substituting in (36) from (37), (38) and (39), we have 


Vv = -V)V+Vicurt | 9) 
2 


V2 . - 


w 
V.Ver 


- wl V+ Vow | =O . (40) 


On examination of the expansion of the last term on the left hand side of (40), it is 
seen that the term produced from V, w is zero, while the term produced from V, 2 is 


by the vector relation 
[A,B] . [C,D]=(A.C) (B. D)- (A.D) (B.C). 


The first term in the curly brackets on the right hand side of (41) is zero by the scalar 
product rule, and this leaves 


| [ Vase | =~ 42) 


Using (42), equation (40) reduces to 
Ww, 2 2 


which is equation (17). 
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Flow about Cones at Very High Speeds 


H. K. ZIENKIEWICZ 


(Department of the Mechanics of Fluids, University of Manchester) 


SUMMARY: Effects of vibrational excitation and dissociation of air on inviscid 
high speed flow past a circular cone, at zero incidence, with an attached shock 
wave, are studied on the assumption of thermal equilibrium. A numerical 
solution of the problem is outlined and an approximate analytic solution for 
the flow between the surface of the cone and the shock wave is developed. 
Two numerical examples are given as an illustration and compared with the 
corresponding solutions assuming constant air properties. 


1. Introduction 


When a body moves through a gas it communicates kinetic energy to the gas 
in its neighbourhood and increases the internal energy of the gas molecules. 
Classical hydrodynamics assumes that this increase in internal energy is instan- 
taneously and equally shared by the various degrees of freedom of the molecules, 
so that the specific heats of the gas are independent of temperature: if it is further 
assumed that poly-atomic molecules have no vibrational degrees of freedom, this 
simple theory is found to be adequate for air, provided that the speed of the body 
is low enough for the temperature of the air affected by it not to exceed about 
800°K; for a blunt body flying at less than 100,000 ft. above sea level this would 
correspond to speeds below about 2,000 m.p.h. 


At speeds above this, the vibrational degrees of freedom of oxygen and nitrogen 
molecules would become significantly excited and absorb progressively larger 
proportions of energy (thus increasing the specific heats) until, at speeds in the 
neighbourhood of 5,000 m.p.h., they would be almost fully excited. 


Further increases in speed would cause physical and chemical changes, such as 
dissociation of oxygen and nitrogen and reactions between them and eventually, at 
speeds above some 10,000 m.p.h., ionisation and emission of radiation; these 
changes would occur more readily at lower pressures, i.e. at higher altitudes. 


The advent of flight at very high speeds has naturally stimulated a good deal of 
interest in the effects of these departures of air from the idealised behaviour. Most 
of the work published so far has been concerned mainly with such elements of 
supersonic flow as shock and expansion waves and with some simple two- 
dimensional problems. The present paper considers the effects of vibrational 
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excitation and dissociation on the inviscid high speed flow of air past a circular 
cone with a nose angle small enough for the shock wave to be attached. Such a 
cone, being a comparatively slender obstacle, would have to travel considerably 
faster to produce in the surrounding air changes comparable to those caused by the 
blunt body considered before. The effects of slenderness and speed are conveniently 
combined in a single variable, the hypersonic similarity parameter, K, which is the 
product of the tangent of the nose semi-angle, 6,, and the flight Mach number, M,. 
For a cone in free flight, the effects of vibrational excitation become significant when 
K exceeds about 4, the temperature behind the shock then being greater than about 
800°K. 


Throughout this work it is assumed that all changes of state take place at 
thermodynamic equilibrium—that is relaxation times are assumed to be zero— 
although a few remarks on the relaxation effects are included at the end of the paper. 


NOTATION 
a= {(0p/dp),}+, speed of sound 
c limiting velocity 
C,,C, specific heats at constant pressure and at constant volume 
d_ relaxation distance 
h_ enthalpy 
K=M, tan 6,, hypersonic similarity parameter 
M Mach number 
m_ effective molecular weight of air 
Static pressure 
q=(u? + resultant velocity 
R_ universal gas constant 
r,@ polar co-ordinates (Fig. 1) 
s entropy 
T absolute temperature 
u,v velocity components in the r- and 6-directions (Fig. 1) 
v=C,/C, 
y’ an effective value of y 
p density 
7 relaxation time 
¢=9 -6,., angular distance from the surface of the cone. 


Suffixes 


c,w,0 denote conditions at the cone surface, just downstream of the shock and 
in the free stream, respectively. 
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Fig. 1. 
Co-ordinates and notation. 


2. Equations of the Problem. Outline of the Numerical Solution 
2.1. GENERAL 


Consider a circular cone (Fig. 1) set at zero incidence in a uniform supersonic 
stream. Then, as is well known, provided that the cone angle is less than a certain 
critical value, an attached shock wave springs from the vertex and, under the 
assumption of zero viscosity and thermodynamic equilibrium, all the flow properties 
depend only on the angular variable 6, there being no characteristic length in the 
problem. Under these conditions the equation of continuity is 


and the momentum equations in the r- and 6-directions are, respectively, 


dv dp_ 


Equation (2) implies that the flow downstream of the shock is irrotational and 
therefore homentropic. 


If we now write 


(2) 


where the suffix “s” indicates that the derivative is taken at constant entropy and 
thermal equilibrium, then a is the speed of sound. Although actual sound waves 
may propagate at somewhat different speeds, depending on their frequency, the 
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definition (4) is consistent with the assumption of zero viscosity and relaxation time 
here adopted. As will be seen later, in the general case a depends on pressure as 
well as temperature; however, if dissociation is absent, the dependence on pressure 
disappears but, owing to the variation of specific heats, a is not proportional to the 
square root of the temperature. 


With the aid of (4), equations (1) and (3) may be combined to yield 


Pu _ut+vcotd 


where v is related to u by (2). 


2.2. EQUATIONS OF ENERGY AND STATE 
Under the assumptions adopted, the energy equation is 
4q?+h=constant throughout the fluid, (6) 


where q?=u? + v?, c is the (fictitious) limiting velocity which would obtain if the flow 
were expanded to zero enthalpy, and the enthalpy A is in general a function of 
temperature and pressure. 


Equation (6), together with the Bernoulli equation 
4q°+ =constant, 
implies that the entropy s is constant throughout the fluid downstream of the shock. 


Since, over ranges of p and T wide enough to cover most aerodynamic 
problems, air can be shown (see, for example, Ref. 1) to satisfy the equation of state 
of a mixture of perfect gases*, pressure temperature and density are related by 


RT 
pom 


(7) 


where R is the universal gas constant and m is the effective molecular weight of air; 
in the presence of dissociation and chemical reactions between the constituents of 
air m is, of course, a function of p and T. 


If the variation of specific heats and changes of composition of air are 
neglected, a? and h both vary as T and are independent of p, so that a’ can be 
eliminated between equations (5) and (6) and equation (5) integrated numerically 
Step-by-step as described by Kopal®. In the more general case considered here 
this simple procedure is no longer possible and the evaluation of a requires the 
determination of the equilibrium composition of ait as a function of p and T. 


*The term “perfect gas” is here used to denote a gas whose molecules are sufficiently widely 
spaced for the intermolecular forces to be negligible; it does not imply constancy of specific 
heats, nor that the gas is chemically inert. 
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2.3. EQUILIBRIUM COMPOSITION AND THERMODYNAMIC PROPERTIES OF AIR 


It is outside the scope of this paper to discuss methods of determining the 
equilibrium composition of a mixture of reacting gases; the application of these 
methods of aerodynamics had recently received a good deal of attention (see, for 
example, Refs. 3 and 4). For the present purpose it will suffice to note that, given 
the knowledge of the thermodynamic properties of the constituents of air, including 
their dissociation and ionisation energies and heats of reaction, the equilibrium 
composition of air at given pressure and temperature can be determined by the 
application of the “law of mass action.” This in turn yields enthalpy and entropy 
of air as functions of p and T and hence p (p, s) or p(h, s) and a’ (T, p) or a’ (h, s); 
these data replace, in effect, the h=C,T, a? oc T and p & p’ relations of the constant 
specific heats theory. Several tables of these functions have been published in the 
last decade“: *: * ”; of these, the most up-to-date set is that published by the United 
States National Bureau of Standards, all the others being based on what is now 
believed to be the wrong value of the dissociation energy of nitrogen. 


2.4. OUTLINE OF NUMERICAL SOLUTION 


With the aid of the tables just mentioned, equation (5) can be integrated 
numerically in a fashion similar to that employed in the constant specific heats 
theory. In a specific problem, the data that are usually given are the free stream 
conditions (p,, 7, and M,) and the cone angle, @.. Unfortunately, it is not possible 
to carry out a solution of equation (5) with mixed boundary conditions—one must 
start either with assigned values of 8 and of the flow variables at the cone surface, 
or with assigned shock angle, 6, and the free-stream conditions. In the case of 
constant specific heats, when the absolute values of p, and T, are immaterial, it is 
more convenient to begin the integration at the cone surface, obtaining a range of 
values of M, and interpolating for the desired M,. However, in the present case, 
when the solution depends also on p, and T,, it is obviously of advantage to begin 
the integration at the shock. 


Having specified M',, p,, T, and 6,, the transition across the shock can be 
calculated as outlined in the next section, yielding, among other quantities, uy, vw. 
a, and sy. With these initial values, equation (5) can be integrated step-by-step as 
follows. The first step gives u and v at the end of the first interval; hence h follows 
from the energy equation (which also holds across the shock) and the speed of 
sound from the table of a(h, s), entropy being constant everywhere downstream of 
the shock and equal to s,. These values of u, v and @ are then used to calculate 
the second step and so on, the integration being continued until v=0, when @=6.. 
The solution should then be repeated for different values of 6,, interpolating finally 
for the desired 4. Having determined the velocity and enthalpy distributions 
between the shock and the surface, the corresponding pressure and temperature 
fields follow from the tables of p (A, s) and T (h, s) or T (p, s). 


2.5. CONDITIONS AT THE SHOCK 


The equations expressing the laws of conservation of mass, momentum and 
energy across an oblique shock are well known and need not be displayed here. 
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Since in the present case enthalpy of air is a tabulated function, these equations 
cannot be solved analytically and it is necessary to resort to a trial and error 
procedure. Eliminating the velocities of flow ahead of and behind the shock from 
the conservation equations, the following relation can be readily obtained 


Pw Pw 
Po 


while the relation between the shock angle and the density and pressure ratios is 


sin 0, ( = - = 9 
Po ( ) 


Pw 


Given p, and p,, it is possible to find from h (p, p) tables the values of h,, and 
pw that would satisfy equation (8) for a trial value of pw. 9%, then follows from 
equation (9) and a,, and sy from tables, while 


Uy = Uy = Qo COS 
(10) 


3. Approximate Solution for Flow between Shock and Cone 


3.1. GENERAL 
In the numerical solution described in the previous section, by far the most 
laborious part is the integration between the shock and the surface of the cone. As 
already remarked in Section 1, changes in the specific heats and the composition of 
air are only of importance when the hypersonic similarity parameter K is greater 
than about 4. But when K is as high as that, the shock wave lies close to the 
surface of the cone and the changes of flow across the shock are far greater than the 
corresponding changes between the downstream side of the shock and the surface. 
Thus, the determination of the transition across the shock, which is a comparatively 
simple calculation, is the more important part of the solution and it seems worth- 
while to seek an approximate method for calculating analytically the small flow 
changes behind the shock. Such an approximate solution can be readily obtained 
by virtue of the proximity of the shock to the surface. Consider again 

equation (5):— 
du u+v cot 


(5) 


| At the surface v=0 and near the surface (9 - 6. << 1), v is much smaller than 
either a or u. Thus, as a first approximation to the flow in the vicinity of the 
| Surface, v can be neglected altogether. 
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Velocity distribution between cone and shock given by exact®) and approximate theories for 
the case of constant specific heats: (a) u-component, (b) v-component. 
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age t 2u=0, 


Then equation (5) reduces to 


whose solution satisfying the boundary conditions on the surface is 


so that v= 2! 
or, to second order in 9, 
u=(1 - 
12 
and v= -26q. (12) 


The accuracy of this approximate solution can be best assessed by comparing 
its results with the exact calculations of Kopal for y=1-405 and 14. Such a com- 
parison is shown in Figs. 2(a) and 2(b) for K > 2:5 and it can be seen that equations 
(11) predict the changes in u and v between the surface and the shock within three 
to four per cent.* Such accuracy is ample for present purposes; however, should an 
even higher accuracy be desired, Maccoll’s series solution™ for u, with which the 
present solution in the form of equation (12) agrees to second order in 9, could 
be used. 


Figure 2 also confirms that the velocity distribution (and therefore the stream- 
line pattern) between the cone and the shock is insensitive to variations in y 
(Maccoll’s series solution shows that for u this is true to the third order in ¢). This 
emphasises still further the governing effect of the transition across the shock on 
the conditions at the surface. 


3.2. PRESSURE AND TEMPERATURE FIELDS 


Having determined the velocity field by the approximate method, the corres- 


_ ponding pressure and temperature fields can be found as already described in 


Section 2. However, a simplification can be achieved at this stage of the 


calculations also, for the results of Kopal (with y=1-405) show that, for K > 4, 


the changes in pressure and temperature between the shock and the cone surface 
amount to less than 7 and 2 per cent of py, and T, respectively (and decrease with 
decreasing y). Hence, it can be assumed without significant error that, in the flow 


_ behind the shock, the specific heats and the composition of air are constant and 


appropriate to the conditions just downstream of the shock. Then it follows from 
the energy equation that temperature is related to velocity by 


T=T.+ 
pw 


*The reason why equations (11) provide a good approximation even when ¢ is as large as 
10° is that at the high values of K considered here the denominator of the right hand side 
of equation (5) varies with ¢ almost as rapidly as the factor 1+(v/u) cot 6 in the numerator. 
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Further, define by a? = . ‘ (14) 
then the presure field is related to the temperature field by 


Yw 


4. Examples 


Because of the dependence of the flow on the free stream temperature and 
pressure as well as Mach number, calculations covering useful ranges of these 
variables for a number of cone angles would represent a major undertaking. The 
numerical results given here are limited to two illustrative examples (calculated 
using the approximate theory of Section 3) which, it is hoped, will serve to indicate 
the magnitude of the effects of variable specific heats and composition of air. 


The cases considered were those of a 30° semi-angle cone at M,=7-860 and 
16°62. TJ, was taken as 218°K and p, as 0-01 atmosphere; on the basis of the 
standard atmosphere these conditions correspond to flight at about 100,000 ft. above 
the sea level (or, in the case of M,=7-:860, where there is no dissociation present 
and p, is not important, to flight anywhere between 35,000 and 100,000 ft.). The 
results are summarised in Table I; the percentage error expresses the difference 
between the constant specific heats (C.S.H.) and the variable specific heats (V.S.H.) 
theories. The quantities given under the C.S.H. heading were obtained using the 
approximate theory and y= 1-405; the values of all the flow quantities at the cone 
surface are within +4 per cent of the corresponding exact values of Ref. 2. 


TABLE I 
M, 16°62 
tit | | 
C.S.H. vam. | | vee. 
T,/T, | 4750 | 4481 | +6 1736 | 12-91 +34 
T/T, | 018 | 1015 | +4 | 1015 | 1-005 +1 
ff, 4837 | 4548 | +64 | 1763 | 12-98 +36 
PylP, | 22°48 2242 | | 9734 | 97-32 +4 
PelPw | 1-060 +4 1054 | 1-030 +24 
Po 23:94 23°76 +1 | 102-2 100-2 +24 
Pyle, | 4733 5002 | 5:606 7-435 —25 
PelPy 1045 +4 | 1-038 1-025 +1 
p/P, | 4949 5-224 —5t | 5821 | 7-617 
34-07° 33-80° +1 33-35° 32:45° +3 
M,, 3-017 3185 | 3352 | 4-200 —20 
M. 2:975 3-149 | 3-315 4-182 —21 
T,, (°K) 10336 , 977 | +59 | 3785 2815 +970 
T, (°K) 1054 991 | +63 3843 | 2829 +1014 
m, | 28-97 28:97. | 0 28:97 | 28-58 +14 
392 The Aeronautical Quarterly 
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At the lower M, there is no dissociation present and the differences between the 
two theories are due to the excitation of the vibrational degrees of freedom of the 
molecules and the consequent increase in the specific heats. As can be seen from 
Table I, this results in decreased temperature and increased density at the surface, 
while the pressure is hardly affected. 


At the higher M, there is a small degree of dissociation downstream of the 
shock (about 14 per cent of oxygen molecules, but hardly any nitrogen molecules, 
are dissociated), but the bulk of the increase in the specific heats is still due to the 
vibrational excitation of the molecules. The effects on the flow quantities down- 
stream of the shock and at the cone are qualitatively similar, though much more 
pronounced. Particularly noteworthy is the large decrease in the temperature at the 
cone surface. Again, pressures are little affected and the pressure downstream of 
the shock is practically the same as in the case of constant specific heats. 


5. Relaxation Effects 


Excitation of the various degrees of freedom of a molecule and its dissociation 
and ionisation are established by collisions with other molecules, each process 
requiring, on the average, a certain number of collisions. Thus, when the state of 
a gas in equilibrium is suddenly disturbed as it is, for example, by the passage of a 
shock wave, a finite time—the relaxation time—must elapse before the distribution 
of internal energy of gas molecules will approach the equilibrium corresponding to 
the new state. Since the translational and rotational degrees of freedom require 
only a few collisions to become excited, the corresponding relaxation times are very 
short and can be safely neglected under the conditions of flow considered here. On 
the other hand, excitation of molecular vibration requires many thousands of 
collisions and vibrational relaxation may be significant in many problems. 


Until recently only crude estimates had been available for the vibrational 
relaxation times of O, and N,. Recent shock tube experiments of Blackman, 
together with the earlier work of Kantrowitz and Huber“, have provided the means 
of estimating fairly accurately these relaxation times at temperatures between about 
1,000°K and 5,000°K. Since the vibrational relaxation time of N. is about 10 times 
that of O, and since almost 80 per cent of molecules in air are nitrogen, the vibra- 
tional relaxation effects in air will.be determined mainly by 7x., the relaxation time 
of N, in collision with their own species. 


From the results of Ref. 9 it is found that 
Tx2=940/p usec at 1,000°K (p in atmospheres) 
= 58/p usec at 3,300°K, 


these temperatures corresponding roughly to the averages between the equilibrium 
and the non-equilibrium conditions behind the shock wave in the two examples 
considered in the previous section. 


Defining now a relaxation distance d=qr, where q is the speed of flow behind 
the shock, it is found that, for the 30° semi-angle cone flying at Mach number M, 
at h feet above the sea level, d has the values given in Table II, 
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TABLE I 
70,000 
= 7:9 | 1:1 25 
d (ft.) = 16-6 0-04 0-21 0°76 


Neglecting the small amount of dissociation present at the higher M,, d provides a 
rough estimate of the width of the shock wave and the analysis (which assumes 
d=0) would be expected to be a good approximation for cones of length much 
greater than d; while for cones much shorter than d, the constant specific heats 
theory (with the vibrational degrees inert, i.e. d=00) should hold. 
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CORRESPONDENCE 


The Design of Intermediate Vertical Stiffeners 
on Web Plates Subjected to Shear 


L. H. MITCHELL and K. C. ROCKEY 


CORRESPONDENCE ON THE PAPER PUBLISHED IN NOVEMBER 1956 


From Dr. L. H. Mitchell (Aeronautical Research Laboratories, Melbourne) 


Dr. Rockey, in his interesting paper on the design of vertical stiffeners for shear 
webs, perpetuates the idea of Timoshenko (Theory of Elastic Stability, McGraw-Hill, 
New York, 1936) that for stiffener rigidities above a certain value the stiffeners remain 
straight and the rigid stiffener buckling load is reached. It is my opinion that this is 
only so when the stiffeners happen to be placed along a straight node of one of the 
harmonics of the buckling pattern of the unstiffened plate. Wiéith unstiffened plates 
under shear loads, the nodes are not straight and consequently the buckling load will 
theoretically always increase with increase of stiffener rigidity, i.e. no buckling mode 
exists such that stiffeners of finite rigidity remain straight. It is suspected that 
Timoshenko was led to his conclusion by over-estimated buckling loads obtained by 
energy methods. 


kL | 


Fig. 1. 


As a simple example of this type of behaviour, consider the simply-supported 
column shown in Fig. 1. It is supported by a spring (spring constant S$) at a distance 
kL from one end. A simple analysis yields the result that the equation which is satisfied 
at buckling is 


P,,t sin t= LS[K(1 — k)t sin t — sin kt sin (1 — k)¢], 
where t= L(P.,/ ED}. 


When k=4 the spring is at the node of the first harmonic and it does remain uncom- 
pressed if S >(4P,,/L)=16x°EI/L*. For this range of S the buckling load remains 
constant. However, for all k other than 0, } and 1, it can be shown that deflection of 
the spring always occurs and that the buckling load increases with spring rigidity. The 
solutions of the equation for kK=4 and 0-45 which are given in Fig. 2 depict this 
behaviour. 
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Dr. Rockey’s Reply. 


I am grateful to Dr. Mitchell for his valued communication. I agree that with 
stiffeners of finite rigidity an increase in buckling load will always be obtained with an 
increase in stiffener rigidity. This has been established theoretically by Stein and 
Fralich*, and more recently by Kleemanf, as illustrated in Figs. 3 and 13 of the 
original paper. 

Referring to the theoretical K/y relationships obtained by Stein and Fralich for 
the case when the stiffeners are placed 0-5d, apart (Fig. 3), it will be noted that for 
values of y greater than 90 the increase in the value of K with increasing values of 
is extremely small, in fact an increase of 600 per cent in the value of y results in 
only a 6 per cent increase in the value of K. 

Experimental results which I obtained for two of the stiffener spacings investigated 
confirmed that K did increase slightly with increasing values of y beyond y,, but 
generally this very small increase could not be detected because of the inevitable 
scatter of the experimental results. It was for this reason that I had to approximate 
and assume that K remained constant for values of y greater than y,. 


Erratum 


August 1957 (Vol. VIII, Part 3) 


In the paper by Professor W. J. Duncan “Analysis of a Vector Field and Some 
Applications to Fluid Motion’ there is a misprint in the last line but two on p. 209, 
where an index “2” is missing. This line should read :— 


“= (u? + v? + + m,? + n,”) (l,u+m,v +n,w)?” 


*Stein, M. and Fralich, R. W. “Critical Shear Stress of Infinitely Long, Simply-Supported 
Plate with Transverse Stiffeners.” N.A.C.A. T.N. 1851, 1949. 

tKleeman, P. W. “The Buckling Strength of Simply-Supported Infinitely Long Plates with 
Transverse Stiffeners.” R. & M. 2971, 1956. 
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